{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np  # 矩阵操作\n",
    "import pandas as pd # SQL数据处理\n",
    "from sklearn.metrics import r2_score  #评价回归预测模型的性能\n",
    "import matplotlib.pyplot as plt   #画图\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('day.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>instant</th>\n",
       "      <th>dteday</th>\n",
       "      <th>season</th>\n",
       "      <th>yr</th>\n",
       "      <th>mnth</th>\n",
       "      <th>holiday</th>\n",
       "      <th>weekday</th>\n",
       "      <th>workingday</th>\n",
       "      <th>weathersit</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "      <th>casual</th>\n",
       "      <th>registered</th>\n",
       "      <th>cnt</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>2011-01-01</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.344167</td>\n",
       "      <td>0.363625</td>\n",
       "      <td>0.805833</td>\n",
       "      <td>0.160446</td>\n",
       "      <td>331</td>\n",
       "      <td>654</td>\n",
       "      <td>985</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>2011-01-02</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.363478</td>\n",
       "      <td>0.353739</td>\n",
       "      <td>0.696087</td>\n",
       "      <td>0.248539</td>\n",
       "      <td>131</td>\n",
       "      <td>670</td>\n",
       "      <td>801</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>2011-01-03</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.196364</td>\n",
       "      <td>0.189405</td>\n",
       "      <td>0.437273</td>\n",
       "      <td>0.248309</td>\n",
       "      <td>120</td>\n",
       "      <td>1229</td>\n",
       "      <td>1349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>2011-01-04</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>0.212122</td>\n",
       "      <td>0.590435</td>\n",
       "      <td>0.160296</td>\n",
       "      <td>108</td>\n",
       "      <td>1454</td>\n",
       "      <td>1562</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>2011-01-05</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.226957</td>\n",
       "      <td>0.229270</td>\n",
       "      <td>0.436957</td>\n",
       "      <td>0.186900</td>\n",
       "      <td>82</td>\n",
       "      <td>1518</td>\n",
       "      <td>1600</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   instant      dteday  season  yr  mnth  holiday  weekday  workingday  \\\n",
       "0        1  2011-01-01       1   0     1        0        6           0   \n",
       "1        2  2011-01-02       1   0     1        0        0           0   \n",
       "2        3  2011-01-03       1   0     1        0        1           1   \n",
       "3        4  2011-01-04       1   0     1        0        2           1   \n",
       "4        5  2011-01-05       1   0     1        0        3           1   \n",
       "\n",
       "   weathersit      temp     atemp       hum  windspeed  casual  registered  \\\n",
       "0           2  0.344167  0.363625  0.805833   0.160446     331         654   \n",
       "1           2  0.363478  0.353739  0.696087   0.248539     131         670   \n",
       "2           1  0.196364  0.189405  0.437273   0.248309     120        1229   \n",
       "3           1  0.200000  0.212122  0.590435   0.160296     108        1454   \n",
       "4           1  0.226957  0.229270  0.436957   0.186900      82        1518   \n",
       "\n",
       "    cnt  \n",
       "0   985  \n",
       "1   801  \n",
       "2  1349  \n",
       "3  1562  \n",
       "4  1600  "
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(731, 16)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>instant</th>\n",
       "      <th>season</th>\n",
       "      <th>yr</th>\n",
       "      <th>mnth</th>\n",
       "      <th>holiday</th>\n",
       "      <th>weekday</th>\n",
       "      <th>workingday</th>\n",
       "      <th>weathersit</th>\n",
       "      <th>temp</th>\n",
       "      <th>atemp</th>\n",
       "      <th>hum</th>\n",
       "      <th>windspeed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.344167</td>\n",
       "      <td>0.363625</td>\n",
       "      <td>0.805833</td>\n",
       "      <td>0.160446</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>0.363478</td>\n",
       "      <td>0.353739</td>\n",
       "      <td>0.696087</td>\n",
       "      <td>0.248539</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.196364</td>\n",
       "      <td>0.189405</td>\n",
       "      <td>0.437273</td>\n",
       "      <td>0.248309</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.200000</td>\n",
       "      <td>0.212122</td>\n",
       "      <td>0.590435</td>\n",
       "      <td>0.160296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>3</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>0.226957</td>\n",
       "      <td>0.229270</td>\n",
       "      <td>0.436957</td>\n",
       "      <td>0.186900</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   instant  season  yr  mnth  holiday  weekday  workingday  weathersit  \\\n",
       "0        1       1   0     1        0        6           0           2   \n",
       "1        2       1   0     1        0        0           0           2   \n",
       "2        3       1   0     1        0        1           1           1   \n",
       "3        4       1   0     1        0        2           1           1   \n",
       "4        5       1   0     1        0        3           1           1   \n",
       "\n",
       "       temp     atemp       hum  windspeed  \n",
       "0  0.344167  0.363625  0.805833   0.160446  \n",
       "1  0.363478  0.353739  0.696087   0.248539  \n",
       "2  0.196364  0.189405  0.437273   0.248309  \n",
       "3  0.200000  0.212122  0.590435   0.160296  \n",
       "4  0.226957  0.229270  0.436957   0.186900  "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y = data['cnt'].values\n",
    "X = data.drop('cnt', axis = 1)\n",
    "X = X.drop('dteday', axis = 1)\n",
    "X = X.drop('casual', axis = 1)\n",
    "X = X.drop('registered', axis = 1)   ###去掉原数据中的cnt，dteday，casual，registered这几列\n",
    "X.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(584, 12)\n",
      "(584,)\n",
      "(147, 12)\n",
      "(147,)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=33, test_size=0.2)\n",
    "print (X_train.shape)\n",
    "print (y_train.shape)\n",
    "print (X_test.shape)\n",
    "print (y_test.shape)    ##讲原数据分为训练和测试"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "# 分别初始化对特征和目标值的标准化器\n",
    "ss_X = StandardScaler()\n",
    "ss_y = StandardScaler()\n",
    "\n",
    "# 分别对训练和测试数据的特征以及目标值进行标准化处理\n",
    "X_train = ss_X.fit_transform(X_train)\n",
    "X_test = ss_X.transform(X_test)\n",
    "\n",
    "#对y做标准化不是必须\n",
    "#对y标准化的好处是不同问题的w差异不太大，同时正则参数的范围也有限\n",
    "y_train = ss_y.fit_transform(y_train.reshape(-1, 1))\n",
    "y_test = ss_y.transform(y_test.reshape(-1, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -3.80696961e-17]\n",
      "[[-1.03249036  0.31165315  1.42065613  0.41928177 -0.03150369  0.06836409\n",
      "   0.04234346 -0.19908101  0.20320023  0.27764841 -0.04960717 -0.08352942]]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "linreg = LinearRegression()\n",
    "linreg.fit(X_train, y_train)\n",
    "y_test_pred_lr = linreg.predict(X_test)\n",
    "y_train_pred_lr = linreg.predict(X_train)\n",
    "print (linreg.intercept_)\n",
    "print (linreg.coef_)          ##得到coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE: 0.182642664743\n",
      "RMSE: 0.427367131098\n"
     ]
    }
   ],
   "source": [
    "y_pred = linreg.predict(X_test)\n",
    "from sklearn import metrics\n",
    "\n",
    "print (\"MSE:\",metrics.mean_squared_error(y_test, y_pred))\n",
    "\n",
    "print (\"RMSE:\",np.sqrt(metrics.mean_squared_error(y_test, y_pred)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAEYCAYAAACju6QJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8FNX5P/DPs5tJsgmQDZeqRAOICIoIEapRalXaiopgFAUpXoqt2n61KmpqVFTwBjVV+63Wr5cqWkUExeanUoq2oFYU5JIAouAFAVlUkJBwyYZsNuf3x2bC7O5cd2d2djfP+/XiBUl2Z87Ohn3mnPOc55AQAowxxhjLLB63G8AYY4wx6ziAM8YYYxmIAzhjjDGWgTiAM8YYYxmIAzhjjDGWgTiAM8YYYxmIAzhjjDGWgTiAM8YYYxmIAzhjjDGWgXLcboAVPXv2FH379nW7GYwxxphjVq9e/YMQopfR4zIqgPft2xerVq1yuxmMMcaYY4hoq5nH8RA6Y4wxloE4gDPGGGMZiAM4Y4wxloE4gDPGGGMZiAM4Y4wxloFcC+BEdBQRLSWiz4hoAxHd6FZbGGOMsUzj5jKyVgC3CCHWEFFXAKuJ6B0hxKcutokxxhjLCK71wIUQ3woh1rT/ex+AzwCUuNUexhhj2amxsdHtJjgiLebAiagvgDIAK1R+dg0RrSKiVbt27Up10xhjjGWob7/9Ftdeey369OmD7777zu3m2M71SmxE1AXAAgA3CSH2xv5cCPE0gKcBYMSIESLFzWOM2aimNoDqxZuwoyGI3n4fKkcPREUZD7wx+wkhcPbZZ2Pjxo343e9+B0mS3G6S7VwN4EQkIRK85wghXnezLYwxZ9XUBnD76+sRDIUBAIGGIG5/fT0AcBBntmhpacHzzz+Pyy+/HD6fD0888QR69+6N/v37u900R7gWwImIADwL4DMhxCNutYMxlhrVizd1BG9ZMBRG9eJNnSKA8+iDc9ra2vDqq6/izjvvxFdffYXCwkJMnjwZp59+uttNc5Sbc+AjAVwOYBQR1bX/Oc/F9jDGHLSjIWjp+9lEHn0INAQhcGj0oaY24HbTMt6SJUtwyimn4NJLL0VhYSEWLVqEX/7yl243KyVc64ELIT4AQG6dnzGWWr39PgRUgnVvvy+h42VSj7azjz44RQiBu+66C99//z1eeOEFTJ48GV6v1+1mpUxaZKEzxrJf5eiB8EnRH64+yYvK0QMtHyvderQ1tQGMnLUE/aoWYuSsJXHt6MyjD3bbtm0brr76anz33XcgIsydOxeff/45rrjiik4VvAEO4IyxFKkoK8HMi4agxO8DASjx+zDzoiEJ9UD1erSpZuZmQmuUIdHRh86ovr4elZWVOPbYY/Hiiy/io48+AgCUlpYiPz/f5da5w/VlZIyxzqOirMSWIWMnerSJDsmbGR6vHD0wKgMfSHz0obMRQuCRRx7B/fffj8bGRlxxxRW49957UVpa6nbTXMcBnDGWcZyYT090iZuZmwn5GJkyZ58OhBAgIhARVq9ejVNPPRWzZs3CiSee6HbT0gYHcMZYxrG7R5tMkpnZmwm7Rh+ynRACixYtwh133IEXX3wRQ4YMwezZs5GXl+d209IOz4EzxjKOnfPpQHJD8nYm53V2K1euxKhRozBmzBjs27cP9fX1AMDBWwP3wBljGcnOHm0yQ/I8PG6Pq666CrNnz0avXr3w2GOP4ZprrkFubq7bzUprHMAZY52e2pA8ADS1tKKmNmAYjO28mcik9e3Jqq+vR3FxMYgIRx99NO666y7ceuut6Natm9tNywgkRObsDzJixAixatUqt5vBWKeU7YGlpjaA6W9sQEMwFPV9n+TF+OElWLpxl+OvPTaZTj5/MtMDscdPh/fwwIEDeOSRR/DQQw/hlVdewZgxY1LehnRGRKuFECOMHsdz4IwxQ+lWOMUJFWUlKMyLH5QMhsKYs3xbSl67k+vb7XgPjQrWGAmFQnjyySdxzDHH4O6778bZZ5+NAQMGWHwlTMYBnDFmKJ0KpzhJK2ktdpzSqdfuZMW2ZN9DO24ARo8ejd/97nfo378/li1bhgULFuDYY4+18jKYAs+BM8YMdZZSoFrJbGoSee1aQ9jy97UmNOVkumSGwJN9DxNdavfRRx9hxIgRkCQJ1113HW688UaMGzcOkQ0pWTK4B84YM9RZSoGqLQnTCjNWX7tWD3ZazfqO76uRl6Ql2wNO9j20egPw2Wef4YILLsBpp52G2bNnAwDGjx+PCy64gIO3TTiAM8YMdZa1zmrryyeXl9ry2rV6sHNXfBP3fZlyfXuyQ+DJvodmbwB27NiBq6++GieccALeffddPPDAA5g8ebKpc9gl2bn6TMFD6IwxQ51prbPakrARfbon/dq1eqphjZVABGBZ1SjD55sdAk/2PTRb/W7ChAn4+OOP8fvf/x7Tpk1Dz549TR3fLsmUxc00HMAZY6Z05lKgdrx2rfl1L5FqEI/t2Wo9v8gnmW5DMq9D6wbg3ON74rHHHsPkyZPRvXt3PP744+jWrRuOPvrohM6TrM609zoHcMZYVkuXtc9aPdjxw0uwYHXAsGdbOXogKl9di1BbdLA/YKLYjNlrYPQ45Q3A66u/wU0PPIZL/v08Whu/xx8Xf47H76tCRdkwaxfGZp0l4RLgOXDGWBZLp/XrWvXb768YYqque0VZCbrkx/e5QmGhOw9u9hpYuVYznnwFvzz/LGxd8EdQXgF+NOFe5JxwDipfXev6fHNnSbgEuBIbYyyLjZy1RHXYucTvi5pfzhT9qhaqLjUjAF/PUq9mZvYaaD1Ofqw8IlC9eBNqn7sTLd9vhv+nl6Pw+DNA5Il6rJvX1ulqdqlgthIbD6EzxrJWtg2nJrLpitlroHdNtm7dgiuueBDFIyeC/L3RY/T18OQVgHLiNxux+9panQLpTAmXHMAZY1krmV3G0lEi+6BrXQN/gWT4uHBwLxo/mo99a94CkQd5/cpQ6O8Nb6Ff93x2STSjvLMkXPIcOGMsa2Xb+nWz+6Ar10EfONgKrye+cMr+5tao+erYa7J3ZQ0CT12NfaveQOHxZ6L31U+h8PgzddvnUTlOMjpLCd9E8Rw4YyyrpUsWeqqozQFriZ2vHjb9X2hojjyv/t9Po7XhW/jPuBK5vfqaOrfkAX7UzWfbtU5kzj8b8Bw4Y4yh8wynytR6rVrk+WohBBYuXIj6l25F+JQrIR11IopH/RrkOTR6IXkIoEjWu5ZQGzqG4QMNQVS+thZA4gVUsm0KxG48hM4YYxlGr1SolSQyf4GEFStW4Mwzz8TYsWPRGmoB2kdllcG7uEBC9SVDUX3xUJRYCJ6hsMCMNzeYfnysbJsCsRv3wBljzAapGqo3SuyysqPaN2/8L8rv+SeKuvdEv3E3IHzsKJA3PiwU5Ea+J7++4gIJe5pCps6xpylkWGhGS2fKKE8Ez4EzxrJaKgJrKtceG63rNpoDDx9ogMfXFeTxYl/dIkgHG1F40gUIS/o9a5/kNT00r/bcmRcNAWBfMM7m3AaeA2eMdXqp2tgi0frbiQQho3Xdar3WppZW7N7TiL0ra7B35T/Q/WfXoMuJv0DXYecCAMyE5USDt/zc6W9swMHWNlvei860YYkengNnjGWtVC1D0gqqgYag5lx1omVetRK4BCK982k166OC99RRR+Pk4CrseOYaNC57Gb5+JyHvyOMTfakJawiGTL0XZrYC5eVlEdwDZ4xlrVRVYtOadyZEZ2Ure4mJ9trVirnIAg1BvLR8W9TXUy67FPu/WIHjTzoF3lMux95ufRN4hc5Rvhdme9bZVmEvUdwDZ4xlrVRtbKGWLU1A3BpmZS8x0SCkLOaipfmbT9B2sAkAUHDSWAy64n58suojrHviOlNZ5LGvxYqR/burZo4XF6hve6p8L8z2rDvThiV6OIAzxrJWqpYhqVVI00oPlgO0VrDxEBkOo1eUlWBZ1SjE1ldr2bUFO1+bge9frsK+2oUAAF/fYWg+YhiIIo9WuyaSl+D3SVHV3awsF5P5fRLmXH2qarW4e8YONnwvzNzU1NQG0NTSGveYzri8jIfQGWNZK5XLkGILxmhli8uBW2soPCyE6YQseei+de8PaPjgJRz4ZAko1wf/Gb9C1+Fj486pPKaZa3LTvDqjlx2lMRjqOIdW2/XOqzUVIc/vnzWoV9ze6UDkxmH6uMGdKoEN4GVkjDHmCDNLy2pqA7hl/lqEVT6HzWzLKZ9j69y7EdxSi64nnY+iUyfA6+sW9bjLyktxf8UQzWNoBdVhM95GQ9Dcem8gEkgL83IQaAjCS4SwECgukCBEJLgb3UBNq1mPOcu3aY5eqE1LAM5tYerWUjVeRsYYYy4y09OtKCvBVI1e7o6GoGYAOXjwIJ544gmMHz8eMy8agnuD12PngVbk+Q9XvRlYunGX6jmMksbOH3qEbkBV8gA40NLaEfDldigLvugt96qpDWDB6oDuubR+Jmf7y9fprEG9sHTjrqQCbyYsVeMeOGOMuUhrqL24QEJzqC2qB5+fQzgn/0vU/O0RbN26FdXV1bj11ls7fm518w+9ojBqQ/wE4LT+3fHpt/tMV2JTo9Zj1mqLGVo9c1kiRXXMFMxxqndutgfOSWyMMeYirUQ7IaKLpwS/rsXmZ36P/73rRnzbnIPjf12NY342Kep5VrOz9ZLG1DLCBYAtu4MdpVUTpXZes0vAYhP3jII3kNgacb1rk+gafrtxAGeMMZuYKUISS2uP78aYueemTR+g7eAB9Bx7Kw6/8lEc6HlcXNCwmnWvF/D1Aliy663VzmtmCZhP8mJyeampbP9YVtusd23SpZAMz4EzxpgGK8OkycyZqmVt3zf3XWx86xl0Gz4WeSWDUHzmFHTPyQPlHFpPHVv4xWrWvdowuRzwqxdv0s2iT3S4W+uGQq0tkpdQmJujmwBndujd6hpxvWujl7eQSq4GcCJ6DsD5AHYKIU5wsy2MMaakFpCnzqvDqq31HRndygDvac+6VjJTWS3WDz/8gPvvvx/r/voEwvDA13cY8koGwZPfRfXxsUHDyv7nRgFfbQ78rEG9MKJPd8tLzGT5kvrAb6JL/vQq08kSWSOu1x6jm5tUcbsH/jyAxwH83eV2MMZYFK054DnLt2FEn+4AogOcWvY3YK1X9pe//AV33XUX9u/fj6uuugrlF1+L52r36vYwY4OG1eQqrYBfUVaCVVvro7LQBYAFqwMY0ae7pS1FlfY0hTRHJqzcfCifA+ivWU90Vzit9uj1zlPJ1QAuhHifiPq62QbGGItVUxvQDJoC6JjrNLNDl1GAvfln/XHh8KPg8Xiwf/9+nHnmmZg5cyaOPz6y4civR2tnlwOIChp6w/iAem9SL+Av3bhLsxzsPWMHo/LVtQi1WV/JZHVkwuimRK9XXOL32b7sK132KXe7B26IiK4BcA0AlJaWutwaxli2k4OgHrO96themTLACiHwxcp38cuHf4Ubbr0D1bf9FlVVVfB44oeYtSqUAYduJvQ2SNHaynPV1vqoymax8/ZGpU1bEwjescfQIgftQEMwKtNcK7cg1b3iREYL7Jb2WehCiKeFECOEECN69erldnMYY1lOLQjG6u33ac53eomissmVH/LysQ8GPsP3L9+GXa/fh7a2Nry0ZhdqagOqwRtQzy6XKZcwaQVFra085674RjebWus1+gsk3P76etMZ4Gr05ouVy7QA/U1hZFrZ/G4HWSelfQ+cMcaUnC5vadQzVPbqjEqlxrY70BBE/b+fxr7Vb8BT6Ef3s/8HXU48G+TNiUuQU1IO2ar1xIOhMG6Zv9ZyQNWat5crm8X2fuXXGLtGPRF6PWMzN1GB9vXYsUPp2RywY3EAZ8xhbtVTzkapKG+pN1xdovL+mXlvZ79Ti1n/2QYgB3lHDYbH1xXdflwBT+6hXqgyQU4rqayirERzPlwrGOvxqmTOA9H7mAscKpYiv36tZVRm+X2S7vtldooi3UqbppqrQ+hENBfARwAGEtF2Ivq1m+1hzG7pUrEpW6SigIbW3t6xzNyY7du3D3fffTd+M+Y07FxeAwAoHDgS/pGTooK3TJkgF3suuUCMh9RaY51P8mLSKUeZ2sdcDt7LqkahoqwkqeVSBGD6uMG6jzF7fDeKp6QTVwO4EGKSEOIIIYQkhDhSCPGsm+1hzG7pUrEpW5jZLzpZyrlUAHEJVFPn1aFv1UJMnVeneWPW0tKCxx9/HP3798d9992H/P4/RsHA00ydP/a1xN4EJtLTlsXOz99fMcTyPuaA/py8HgIwubzU1Npus8dPdfGUdMJD6Iw5KBUBpzPRGt62u4CGPFytVuVLxPwtUy6NmjJlCl5++WWceeaZeOihh3DTf/aZrlxW5JOivjYzH2wGAXh4wlDDtddG+5jLz5HbptwBbN7KbxAKR18Zn+RBc6jN0vSR2vEPHGxV3do01cVT0gkHcMYssDqfbUfA4Tn0Q1KxVEh5va30dZu3rsO2/UcBAKZOnYrJkyfj3HPPBRGhMid+b3CvhxBWWYbVEAxhWs36jmQ2u272BMzNFZu5xmq/k9WLN8UFbwBoDrXh0YnDEirQonyO1v7qqS6ekk7SfhkZY+kikflsq5tL2HHObOb0UqHY621Gy64t+P7Ve/D9K3egbf1bAIARI0bgvPPOA7XPV6u1++FLhqK4QFI95pzl2zreY7t6mCUmj2N0jdV+J+XpBDVa8/pWdcZlYkZ4P3DGTDLaH1hLMj3oRM/JEmNlT+rWvTvR8N85OPDJEnjyCtDjJ5fisfuqMPHUY0yfT6/CmvweT6tZj5eWbzN9TDWJ7IetJZF9u7X2I2fqzO4HzkPojJmU6Hx2MmtTO8McejpNEVi5rg3/fQkHPvsvup18IQaOvhx3XPhjAJEAl+wUi7ItSzfusvAK1I0fHqnSNnVeXdLXOJHfvc48T+0kHkJnzCS9/YHtJi8b0uqdZcsHYrpNEehdV9HagsYVr6Nl59cAAP9Pr8CA6/6GF578Cz6+twIAVIeWp9Vol2WtHD1QdYkaEBl6Nurtys8t8fs0h+MLJA/mLN8Wd42n1ay3vHc5YP13r7PPUzuJAzhjJiU7n21WbBnJWNn0gZhuy+zU3mPRFsb+9f9B4Olr0fDuc2j6/EMAQE7Xnmjx9ehYVnbL/LWau5dpBceKshJMLi/VDOJGQ9WPThyGLbPGYFnVKNwzdnBc2yUPoSnUppox/5JKUNcL4vJNpVydTY9eOVlmHx5CZ8ykVO1ApLdsSK0SmBNSNaxt5xRBsm2Wnx8MhTsqlB3cUovdS55FaNcW5B4+AD3HTEV+nxOjnicHR6312XISl1Zb5GzzROa5Y8uIAtG/n00traa3/FQug4u9lmcN6hW16YmyOptaqVUO2qnBAZwxC1JRa1kreBGQksS1VJQrlQOE1hSBv0CyNJecbJtjnx8WAj7Ji5LWrdgdakbPcX9AwaCfgCixQUujG5JE5rnVsspjfz/7VS20dMwd7fXFY6+lck9wmbK0qpkbJydvCtMpjyKVOIAzlmZSVaxEi96wth0fimrreZUkL2F/86Geo5lgPP2NDUm1WX5+aM+3aHj/7ygcfBZwzMkIn3gh5vzhDjyw6AvVIiJmGb13VkccJA+hoakFfdsDtN8nYfq4wXGvVS9JTqudau+/XnU2Mze1Tt4UpuKGM13xHDhjaSZVc+1anM58N5oiKMzNQSimwInevHhNbUAzuMq7ahnN7e7+YRfq33kSO/72WwS/+hjhfT8AAL5vErjk5H6ou+ds/HniMPh96oliesy8d2ZuzuR5Z79PQhuAAy2HrmFDMITKV9fGvU4rJUnldlp5n83eVDqZ65BueRSpxAGcsTTjdsEKp7Pt9QLEjoagZjDWep7RB7VRgtbN0/+IwNNXY1/tP9FlyC/Q+5pn0LXsPADx5UPlQG62KIqXyNR7Vzl6ICSvfmqYPGRdmJejWsEt1CY098hWu/GQvAS/T4r7HdN6n2NbZ+Wm0smbws6w1FILD6Ezlobc3Nf4rEG94uY8fZIXZw3qpTsvbXYeUm9YV6+slFpgmVajna2vpOyRVS/ehED9fhxRlI/bzhuMvQcF8vsMRfEZV0LqcVTU85Id9WgTwlLt7xlvbtBNOjMKSmo/l3+XzL4/WqVUxw8vwdKNuxKaZ3ZyWsjtKSc3cSU2xjo55Qd7kU/C3uYQYjt4I/t3x5ptjXEf6nKvTatOtVrvs6Y2gKnz6izVGY/dj7qirCShCmX5OR7Uf/oB9rz3d3Q9aQwOK69ArpfQ2Nwa99jiAgm1d58d13a9+ftYiVbM06vAp5dZbleFPruTwqz8fqTTsd3CldgYY5rkD2h5Ta8cTLWGr5d9VR/3PWWSmJXEt4qyEtw0r850W2O385QTlOau+Mb0MQCgefun+O7d2TgY+Aw53Y+E5D8CwVAY+ZIEn+SNCwD3jI3fs9rKzmCJ5i3U1AZw4GD8DYV8vOlvbNB8rl15EnaPADm5BDNVyzvTEQdwxjqZ2B5LMmNwci9Ra2hXTiJTfqDW1AY61lnrkeeZY3ui8o2BlX2x97z3PPYufw3eLt3RffT16HLiL0CeSHLXnqZQexUzgWCoDQCQL0WnB0Wu2bqOn6vx+yQQAQ1NoaggYqU3q9XDLy6QcM/YSJb5VJ2bn3QOWk5OC7k55eQmwwBORIcBeBBAbyHEuUR0PIBThRDPOt46xjoRp9ayxh63qaXVlv2lAYAIuiVfgUMlRW+aV4fiAgn7m1sNg6/c29QKVoGGIDyEuKF+pdZ9u0FSHrz5XeDrexI8kg9dR1wAT25+3GNjh6T3NIU6evoAcPO8OmiHbu2ha6tLnLR6+AW5OR2P15rzNZtYx7KH4Rw4ES0CMBvAnUKIoUSUA6BWCDEkFQ1U4jlwlq2cmsezOmebDpQ9WUB7hMADqAbVk3vnYcn8Z1G/4h/oWnYeikf9OuG2aI0CKEleQmFuDhqDobgbL625bC8R2oSIe7zW7mTK3byycc6XRTM7B25mGVlPIcR8tP9fEUK0AsicTwPGMoBTa1mtzNmmg+ICCQdb27CnKQQB/eH9NgA+yQNv+57bnnArjt31PpbeNwm7l72CHoNHostJyW1huaMhaJj5HQoLNARDqjXFtZ4bFkL18WaW8Lm9zJClDzMB/AAR9UD7/yUiKgfQ6GirGOtknFrLavb5yl2tLisv7QiKqbanKWTphiMYasPhRfkgAHvffQbvPPcQ9heW4MTr/w/PzP47cv2HJ9We3n6f5eVIyhsvM89VPt5sEZ+KshIsqxqFr9s3MuHg3TmZSWK7GcAbAPoT0TIAvQBc7GirGOtkrK5ltWPNtdKjE4dFPX9OAhtruKF5Sx22FP0IUnFv5JWNw4/6nYz8fiehkQi3v74e/gLJ9GYesZSB02gOPJZ846S2plqN/B515oxqZp2pdeDt894DEblR3ySESLwocBJ4DpxlK6vrqNWCglo9bDNz4MoELOXyslST55LN1Bxv2bkZe959Hs1fr0GXoeegxznXqz4udqcsLV4ilB9djC27g6qB00wWupLaNd3RvnWn3nM4WDPAxnXgRHQdgDlCiA3tXxcT0SQhxBM2tJOxTiu2F2220pXWvHZDMBSX4azs0cWu+Qaie5luJ7wV5uZg+rjBcW2QPIQu+TnY0xRCa+NONPz3RRzY8C48+YUoPusqdD3pfM1jml1oFhYCa7Y1as4lK5cpaSWmyWKHvJXP7auzO1hn2oSD2cNMFnqdEGJYzPdqhRBljrZMBffAWbZIJpNYK1NZpleNK7bqmnLdspW9o51CQFy75BuZkbOWYP3rj2F/7T/Rdfg4FJVfDE9+F1vP7yXCwxOGGm5dqlVJzuj5RsEfsK+amhmddRvOdGdnFrqH6FBGCxF5AeQm0zjG0kVNbQAjZy1Bv6qFhrtW2SmZrPNktqaUk58enTgsKts70BB0PXgDkR5zQzCE5lAbHp04DP++8VR8/s4cvPfee5GbjlMnovfVT6H4zF/ZHryBSE9cb+MTIHINtW6gjGqfm9kdLFWbcMg3kYH2oX2jTV9Y+jETwBcDmE9EPyOiUQDmAviXs81izHlufoDpVS4zupkwCgJmMp/TfXlZ08EW/ObOahQdXorbbrsNv7z7CRT5JHh9XZHTrZej5zZzI6VVNKXIJ+neECqXgGnxEKXkd7Azb8OZLcwE8NsALAHwOwDXAfgPgD842SjGUsHNDzC9IGt0MyEHgUj5z2hm62+n81aLwa/X4NvZv8fuRf8LKuyOwybNhHTqFTjQEl8f3ClG10ftJkryEA60tBreEMqjIH+eOEz1RszMKIAdOvM2nNnCMIALIdqEEP8nhLhYCDFeCPGUECJ9b90ZM8nNDzAzQ6lqNxPykP/UeXUoyM3BZeWlpgt6KKcLPCle522mpyAL7f4GIhxCzwuqcPjlDyO/NFL0MRRO3c6JRqMYasVUuuTnxLVR74ZQPobamvtU3Eg6ve87c55mFjoRzRdCTCCi9VBJ5hRCnOhoy1inlorkGjf3EY5d76sVmpQ3E2p1tResDphKfIt9rpWNQLRIXjIdVIsKJBTk5qhe71B9AA3v/x2+o0egy4m/QNeyMehaNgbktW+vpeICCc2hNlPTBmZHMWI30OinkWFulJOgVe/d6RtJrX2/7drRjDlP78b4xva/zwcwVuUPY44wmpu2K/HMbNUrpyiraWnNiSpvJpIZ8tea87bSEc/1Ukdv0UuEiT8+yvQGGg1NoY5hY1n4wB7sfvv/sOPZ/0Fw82q0hZojbfLm2Bq85a1Bxw83vgFMpixpoj1at3rCXJI182n+LxFCfNuecf6sEOLnKWwT6+SMApWV3Z30pKrqlZnRBDO9Ia3lR/KWnXrH1+zNCeDPE4dh+hsbDAuotCh622EhMGf5NgiYK5YiB6OKshL84bW12L16EfYsfRYidBBdhp0D/2mT4O1SbHAUdWrnl78nF0cBgAWr9W/0CEhq+VaiPVo3e8KddRvObKF7myuECBNRExEVCSG4/jlLCb25ab3gnsgHkdMfYGa3kzRzM6G3h7Yc3LWOrzddUFFWgurFm0xVQFMSir/lgOn3STjQ0ho1tO6TvDhrUC+c+sDbCNTvh0fKg7fQD19wbCKsAAAgAElEQVTfMvjPuBJSd+vXv8Tvi7pOgP61GzlrieHwebI93kRvCLl8KkuUmXGqZgDriegdAAfkbwohbnCsVaxT0ws2mZY5a+WGw+hmwuy8tdrxjXp5yV4/ube7rGpU3IjDmQN7YvZL87Bz6fMoHPRT+E+fjIIB5SgYUJ7w+ZpiMtKNrp3R65O8ZEuPN9EbQu4Js0SYCeAL2/8wlhJ6wUarTrfbmbNaw+R23nCUmNyYRO34er28mtoAPDq9e7N2NASjaql7ifDV+pVY9dhstOzYBKlHKfJ62zMsLBedMTuForepS3GBhHvGDuYA6iCu+OYMwwAuhHiBiHIBDELkRnuTEKLF8ZaxTstoSDHdMmf1hsntzHSvHD0Qla+tNZX5rXZ8tV6e3HY7stJ9kifqOuz+4GU0fjAH3i490OPcG1B4ws9AHv2lc4kwM4WidVPISVvOMzuNxKwzs5nJeQCeAvAVIlNd/YjoWiHEIqcbxzovrSHFdJwv1BsmVwu6SQ3XxsRZDwBvzHIuoxsaZW/Ijp63rCnUhta9PwAeD3K6dEfBMaeAvDnoOnwsPFK+LefQYpTIl46/N52F3Xkr7BAzQ+iPADhLCPElABBRf0SG1DmAM1ek23yh4TB5bHxMMF5WL96EUFv0k9sAdMvNQWFejqnAVFMbiLqhsBq85Xnusnvfjqqd3ta8H40rXsO+VW+gYNBP0XPMTcg97GjkHna0tRfZTvIScjxkevtOgnEiX7r93nQWmZa3kknMBPCdcvButxnATofaw1jG0RsmVwu6oTaRUO9D6wOvMRhC3T1nmzrGjDc3mBqC9xAQ0+yonn1De/AWrSHsW/MWGj+aj7bmfSg8/kz4fzLJVFv0suon/vgo3F8xBIPv/hcOtBgXX4k9SmwPr6Y2ELVUjue9U8fNgknZzkwA30BE/wQwH5H/J5cAWElEFwGAEOJ1B9vHWFpRS8bRS7qzs8qWHR+EZncciw/eno75YmXSW8MHL2HvigXI71uG4jN/hdzD+ptuS5tO73/uim861pknSr7GNbUBVL66NupGak9TCJWvrQXA87BO44pvzjFTojgfwPcAzgBwJoBdALojUo3t/GROTkTnENEmIvqSiKqSORZjTtOqEAdAs6KVnVW2tCrHnTWol+Nboja3D2XX1AZwQ/XzCH4XGZTrOuIC/GjCfThs4n2WgjdB/xqEhUgqeENxfLVRECBSW5133nIeV3xzjpks9ClOnLi9yttfAfwCwHZEevVvCCE+deJ8jCVLLxlnWdUo1Q8kO3sfaolYfXv4onqqypuK2MdWjh4Iv0+yXLAFiAy93TP7LWxd9DQav1yDguN+il7j/oCcLt2R06W75eNNLi/FiD7dMXVeXdKBWo3ZNe48D5sanH/gDPsKDlt3MoAvhRCbAYCIXgFwAQAO4J1IJq0PTSQZx+7sZ+UHYU1tQDUABkNhTH9jAw62tsUt3TmptAjLvqq3dM5Qw3do+O+L2Prpe/D4uqF41NXoWnZeQu33EmHSKZH57ZraAPIlj+lENZkHkc1RtKYDSmKusd4acDfmYTPpd56lNzcDeAmAbxRfbwdwSuyDiOgaANcAQGlpaWpaxlIi09aHJjoH7VTvo3rxJs3eq1ovOxgK40OLwRsAmj57H8HPl6Oo/BJ0K78YnrxCy8cADmWwA8r3Pjp4qyXPxSpqT0Azu667cvTAuDlwwL7qa1Zk2u88S29Wtum1m9o+SGrblj4thBghhBjRq1evFDSLpUoyu2u5Idndy+zaRU2WyPCv0XC1T/LilCML0PjRfDR9/iEAoOuIcSi55in4z7gy4eAdGyy1dkczs6qtoSlkaV61oqwE1ZcMhd8ndXyvuEBC9cVDUx40M+13nqU3vf3Ab9Z7ohDikSTPvR3AUYqvjwSwI8ljsgySaetDkxkOt6PnFTv0WpTgfLYW0RbGKPoUr9z7MBp2fY+uJ52PgmNPQ4+irthjsRCLcq5dbcmW1ntsZj5cubOZ2WuXLnOwmfY7z9Kb3hB61/a/BwL4MYA32r8eC+B9G869EsAAIuoHIADgUgC/tOG4LENk4vrQRAOBUTUqo3lRtRsAIDKEZm0GWV1wSx32LnkaT+zahtzeA3HYL6ci/6gTAEQy0It15pzVNDbrP1ZvXlopdqtQN5Yf2TlnnYm/8yx9aQ6hCyFmCCFmAOgJ4CQhxC1CiFsADEekt5wUIUQrgOsBLAbwGYD5QogNyR6XZQ61IWnJQ2hqaXV0SZQb9HpeWsvTlK9da8g52eAt2sesqbkR4XAYvSruwOGX/akjeAORGw0hEPde6R/30L/lNdfK16P23qseB3B1+ZGZ98aKZKdhGFMyk8RWCkC5eUkLgL52nFwI8U8A/7TjWCzzxA5JF7XvJW11pym7JNLTMvsco2ptRrWi7R5iDe3ejob3/47cI47F8edcjv1DRyFvwEiQV/0joTEYwqMTh2HGmxss9cQ7zheOrj4X+95r1WRXJr65wYn95+XjchY6S5aZAP4igI+J6B+I3BBfCODvjraKdRrKIemRs5bEzemmatODROaorTwn2WptZoecjYT370HDspexf+1ikJSHvJLjOtqgFbzl81eUlaB68aaEAjgQvd2o2pam6Vity4k563SZj2eZz0whlweIaBGA09u/NUUIUetss1hn5GaCTyI9LSvP0ep5AdDsfSrnRStHD8RNGoHerP3r/4P6d/4PIhxC17LzUHTapfAW+lG9eJPhDYKZoihGinxS1EYqgYZgXDnTdOuZ8pw1S2dm14EXANgrhJhNRL2IqJ8Q4msnG+Y0LqaQftz8sEzk5sHMc/R+z/T24o7tfa7aan39NgCIcAiiNQRPXgGk7r3h6/9j+H96OaTi3lHtfXTiMMMbhJGzliRcNc3rITQGQ3HPD4UFZry5oaNX6tb/Qa33iet4s3RmZj/wewCMQCQbfTYACcBLAEY62zTncDGF9GTnh6XVG7REbh6MnqP3ewYAt8xfqxq8vURxyVpzV3wT9zg9Qgg0bfwADe//Hfl9h6LH6OuRV3IcepUcp9reirISzfnt4gIp7n0xosweL8z1oqW1DVrP3tMUQt+qha7tEGbm84Bv9lk6MtMDvxBAGYA1ACCE2EFEXfWfkt54g/n0ZNeHZSI3aGrVuiSPfqUuoxsOrd+zqfPrdAuWtAkR104r+3Y3b1uHPe/ORsu3X0Dq2QcFx5RrPlbZXq3qZkLAUvAGIjcFcvLZyFlLTM3fy9nqq7bWY+nGXSkLmEafBzxnzdKVmQDeIoQQRCQAgIgSK8WURriYQvqy48My4Ru02NqAarUCFeRjKXuueTmHVmZqFisxiMVqvX69vbOV9n78OvYsfQ7err3Q47ypKBx8JsijvVxL2dPXuoHSSrLTo3ztVv5fhcJCc3MWp4Iofx6wTGUmgM8noqcA+InoagBXAfibs81yFiemZLdEPpCrF2/qSK6SyfOzRiMCzYp63g3BEKbOq8NN8+rgIXOlQZW0pgwmnXIUXlq+TfU5rXt3QYRbIRUfgYKBIyGEQNeTzodHytM9V3GBhOrFmzB1Xl3Ua4t9fdWLN1nOgFf+X7KaQa+2OYuTo2P8ecAylWEtdCHEnwC8BmABIvPgdwsh/uJ0w5zExRQyi9Ua4onswa0V3Pc0haKKeFS+uhZl977d0Zbpb2yI6+3LAchoU45YanPfshF9usPriR4SCDfvx56lzyHw9DXYs+QZAEBO0WEoOmV8VPCWvAQp5rmSl7C/udVUgRKzRVdin5PM82M52RvmzwOWqcwksf1RCHEbgHdUvpeRODElcyQ6n201Gc5sLzHUJqIKzdhFayctWfXiTQi33xGI1hbsXf0W9i6fj7bmAyg84Sz4f3KZ5rGrLx7acQz59/3AwVbdNfexSYDjh5dg6cZdCDQE48qbqlG+DuX/N7PPj+Vkb5g/D1imImEwxkdEa4QQJ8V8b50Q4kRHW6ZixIgRYtWqVak+LXORVgKUUYUuq1noaoVEUsVM9nW/qoUdQa/xo/mR7PJ+w1F85pXI/dHRms+T98aOvRZ6S8b+PHFY3LUgAJPLSzv28dYbVtd7b/QS2ooLJIw58QgsWB0wtU0oY9mKiFYLIUYYPk4rgBPR7wD8D4D+AL5U/KgrgA+FEJPtaKgVHMA7H2XgUiIAX88aY+u5lIHJbNKYHbYYvA4hBAZf9UfUHxTILz0RbQebcPC7L+DrM9TU8WN7vHo9YC8RDi/K1wyyl7UHcSByvZSFWWSFuV40tYRVb5zMvJ9co4F1dmYDuN4Q+ssAFgGYCaBK8f19QojEqkowZpGVBKNkP/jlx6ayJ16Yqz83/Kc5/8R9d9+JvZvr4Ov/Y+SXnghPXgF8fYZC8hIgELX0TU3sT/UeHRZCd755zvJtGNGne8e1Cofjj3agRXu6w8z7ycu2GDNHbzeyRiHEFgD/C6BeCLFVCLEVQIiITklVA1nnZjbByK5do7R2/XJKU4v6uTZv3ozTz6lA5WVjsP/bzSj++bXodeEdHSvbSvw+VF88FNWXDAUZLHezosTv051vFohcI7T/bbQbmjyvLuOEMcbsY5iFDuD/AOxXfH2g/XuMOa6irAQzLxpiuKWk3tpvK6xkOxcXSPD7JEvHj6XVG16yZAk+XLIYRadORMm1f0O34WNBXqlje81lVaM6eqp2jfTLgfSsQb10HydfI7PXSvk4s+8nY8yYmXXgJBQT5UKINiIyW0OdOSyV84VuzU2aGVJNtBhH7Gsq8klx2dkA4PdJKMzLMbWLViyf5EVza1g10BIiiV3bd+6B+GQhxpYfh7/eW4niob/AEdc8jZwu3XVfU01tIKGs7lglitdkdNMj99DNZu57iNCvaqHuWnPGmHVmAvFmIroBh3rd/wNgs3NNYmalsqZ7sudyOvgnUoxD7TXJa6aV88o+yYvp49SzxCvKSrBqaz3mrvgGYSHgJUL50cXYsjvY8VrPGtQLr6/ejqZQ/IBzW1sYG9/9BxqXvYzw/nrM2XUOik86BwtWB1SDd+xrql68yVTwloN8bLBXy/DWu+lRDnerlZ9VIycD8p4DjNnLzBD6bwGcBiAAYDuAUwBc42SjmDl2DRs7fS675qf1WJ1brakN4Jb5a+NeUygs0CU/x/QQb01tAAtWBzqCVFgIrNnWiMrRA/H1rDGoHD0QC1YHVIN389Z12PHsdahf/Dhyig7DYZMfgn/09Zi74hvNHj0hcv3kgjZ6wVb5Gh6dOAxbZo3BoxOHGb42vZse+T2vqQ1EeuuXDI2aRiiQPCgukECIZLRrPZ8xljwz+4HvBHBpCtrCLEplDedkzpWKzWOsFOPQ28YTABqaQqi9+2xT5zV6bWo/F21hkMcLISJBvdeFd8I3oBzUHvD0lq/F1gj3F0iqO4hprcU2M3ytVghHKbYnrXW8flULVb8v/86YGZXhJWWMadMM4ET0ByHEQ0T0GFSm2IQQNzjaMmYolTWckzlXqm40zM6tGmWaq83ZajF6bcqfh374BnvefwE5/sPRfdRv4Os7DL1//de4zUbIZA31YCiMvBwPfJLX9i1Yi3wS8iUPGppC8KisiQ+Gwrhl/tq4OupKer8zZqZkeNtfxvTpDaF/1v73KgCrVf4wl6VySU4y50qkNrkaqzXRtRjdOISF6BjqnzqvDn11zmf02nr7fWjdtxu7//UYdjx3HZq3roW3sLjjcbHB2+sho03QojQGQ7ZkdcdOczQEQ2gOteHRicPQpnE3obxOalMier8zZqZkUjlFxFgm0uyBCyHebP/7hdQ1h1mRyhrOyZwrkdrksezsjVnZHctoW0uj13bYD6vx0TN3QYTD6HrSGBSddim8BUWa52trE+qVyjR65b39Pse3YDVzvYKhMGa8uUGzBnrs74zWFqVmtiHlbT4Zi9AbQn8TOqtThBDjHGkRs8TKh7cdlcoSCRR23GgYzTVbeW1aQdeogIvavL3aa7tpVD+MPDIXNbUBrN5fBN+AcvhPvxyS/3DD16n1H06I+DaauQkye130guWjKrXR1expCnUkt8m0fmfMTMnwNp+M6dNLYvtT+98XATgcwEvtX08CsMXBNjEHODGfaCVomg3+WsfUCzBWX5vWDYWZfa/l8ynniokiiW9HdMtDOTaiatLVGDp0KL4/+XqQvwS9xlbGHUer1rrW97U2JbGyQYveddELlrHXS21OXGY2MdHMqIwdIzeMZTMzu5G9L4T4qdH3UoE3M0lcort6aVErYJLsrlFqx5Q8hC75OaqZ1kAk4HXzqf/c6mszU5SluEBCc6gtfjRg61o0vDsbLd99ib7HHo8n//IIfru0VXeNtlqP+qTSInz4Vb3hWm0zrLznVt7PmtqA5m5mVjaZ4Sx0xtTZsZmJrBcRHS2E2Nx+4H4A9GstsrRj93yiE0vD1I6p3H9bTVho/1xeL232w99o32qf5IUQiGvjvrpFqF/8V3i79UKPMTfjiNPOxejRP0fvWu2tM+Ue9Yw3NyjaL7D86z1xQX/8cPXRC6PgZuU9l5+nbE9ejnqOa0VZCaa/sUG1Yp2V4W0zozJctY0xbWYKuUwF8C4RvUtE7wJYCuAmR1vFbGdXJrhMKzAlk2Bkd3KSXPTESvGYirISLKsapVn0pLE9aLXu3YmWnV8DAAoGjkTxqN+g5Oqn0OWEUfh270EAkSFgyRufUy55qGMYuFlR4CUYakNYparZwnXfxn3PTHEcK+95TW0A09/YEHUz1BAMaV6z6eMG86YkjLnMMIALIf4FYACAG9v/DBRCLHa6Ycxedi45k+tvq0kmwcjO5CS1+uBmliDV1AZQdu/b6Fu1EDfNq8OBg614dOKwjs1DeuWGsGfJswg8fS3q345UF/b6uqHbjytAObkAIuvI5aBXmBs9yOX3Sai+ZKhmkRc1aiMMZpZYWd3JTa1HrXXNeFMSxtxnOIRORAUAbgbQRwhxNRENIKKBQoi3nG8es4udS8606m8TkFQPzKgCmB55s5FAQ1AzEQzQ7+VPq1mPl5Zvi/peQzCEylfX4mBzEFv/+zo+/fODOLBvLwpP+Bn8p09WPU5YCFS+uhagSGlWGbUfTw6IVkYcYofLtUZAlN83+54b3UhotZOHtxlzl5kktnmIFG65QghxAhH5AHwkhBiWigYqcRKbu+QgopepvcVkApPROXY0BJEveRCMqSEueQkQiNtsZOZFQwDA8AZAK7GtpjaAqfPqNJPOvJv+jc01f8Z5552HkrOvwjvf5ie1A5jkIRTm5aj2emMVSB4IUNTr0tqBjAA8OnGYpWSwflULdV9LoomOjLHEmE1iMzMH3l8I8RCAEAAIIYKApWJRLAso51y1lNgwBC7PQT86cRjUfs369ihAl/xDA0d+n9QxdGvUk9SbMogdVRBCoOmrlWj6YjkAoPWYM3D4L2eh4fRbcN5Py6PmxxMRahNoaQ3HDXHHkjyE3Jz4Neqa68XbX4uS0Xy53tQFz2szlr7MBPCW9l63AAAi6g/goKOtYmknmeBo5/m+2Hkgak74YGsbVm2t11wyJTOao1UOEx/csQnfz70du16bgb0rawAA5JWQd9QJCDQEcdO8Osx4c0PHjmOJ3rg0hdri5pEvKy+N+rr6kqEdiXNmxQ55G82Xq82VA5ElczyvzVj6MrOM7B4A/wJwFBHNATASwK+cbBTT5ta6WKNtK+1uh9n54WAojDnLtyU9BNzb78OWrzej4b0X0LTpA3gKitD9F79Fl6GjVR+/pynUURQlmbl7M/PIWku2tMT2qI2Wk6WyJC9jzD66AZwi+xtuRKQaWzkiY5o3CiF+SEHbWAw3d2fSSpxyan40kXrlaswm1lWOHohr738Xwc2rUHTaJHQ7+UJ48gp0nyPvyPXwhKGYedGQqADY1NKqu34dADyEuNKjqq/B4jh97OvVupZFin28OSGNscxjJolttRBieIrao6uzJ7HZXU3NCicqrxmdTy+pzAplYp1yBOMwH9Dvu6UY2qcnqqqq0Oe2t9AW3Ku72YgateugV63M6LmxjJLMlPw+CXX3RO9lXlMbQOWra6MS/5ScGEFhjCXOziS25UT0YxvaxJLk5u5MqV73W1FWgsnlpYZJYkY/V85Pyzch23fvw97af2Jl9WV45cmHsXhZZHfcI4sLLAdvQH2tdEVZCYoLJI1n6D83llaSWexr90leTB83OO5xFWUlUYl/sdSKwNi1dStjzDlmAvhZiATxr4hoHRGtJ6J1TjeMxUu2mloyH8puzL3fXzEEj04cBq/GGLKXCJPLSzUzuWMT66oXb8KezWux47nrUP/2E5CKS3D4ZX9Cy8jfAdCunOZBZCmXHrWbqHvGxlcrM/tcJa2CLJNjEt70bqgaDIbzlTcSZqq8McbcZyaJ7VzHW8FUxQbNswb1woLVgYS2lIxNhLIyfz6tZn1Uolgq597l48cO3xMiBVOWbtyF8cNLsHTjrqgiLrHDwqFQKLKTVm4ByJODXuPvgq//ySCiuGQuZT1wv0/C9HGDO7YsvWX+WtUiMWo3UWZ38TK6AbMjycxMToF8HZyoc88Ys5/mHDgR5QP4LYBjAKwH8KwQojWFbYvTmebAteac5WCV6JaSsYzmz/XmolNZ4EOviIzkJVRfPFT1Onz66aeoqqqC3+/HV8f/KtKrFAKk6NVbeR3J5AKkOo/A6Nyx5OugNeduZacxxlji7JgDfwHACESC97kAHrapbcwErV7Q0o27sKxqFL6eNaajPrfV4ygZDd9qlU0181w7yQVe1OaVQ2GBGW9uiPre9u3b8etf/xpDhgzBe++9h8GDB+PWs4+FT/JGBW+r69eTyQVws3648tyA+vy5fB3s3viGMeYMvSH044UQQwCAiJ4F8LFdJyWiSwBMB3AcgJOFEJ2jW22BXQlrRo/vHZPkFTtMq/f8RD7Qk51L11qatacphJGzlkSC0NaVmDRpEsLhMG644Qbceeed6NmzJwCAiJKey7ey5Ert9bpVllTZbr33QW1dO1dkYyz96AXwjk9KIUQrWV2Mqu8TRNaWP2XnQbOJ1pyl1aCpN/ep/FDWWmNe5JNUi4gksnGJk+vYRWsIW7f9gNtfb0Hl6aW49NJLcffdd6Nfv35Rj0vlemc31+0b0bsOXNiFscygF8CHEtHe9n8TAF/71wRACCG6JXpSIcRnAGDzTUFWsasXpFUlrLhAwj1jB0d9WKsN2edLHvgkb1wC2eTyUssf6HYkR/ljbiiEaMOBT99Dw39fQk7XHvD+8o94dk0Dls2ebaltTnA7GSyZ0Q4u7MJY+tMM4EII4/UvKUBE1wC4BgBKS0tdbk3q2NULMnscraHyhqYQHp04zJbemNY5Ag3BjuFvw7Ki4wZ3FCUJbqlDw7uz0fL9V5B+dDSKTpsUlVXuNjfX7adz758xZg8zy8gSQkT/BnC4yo/uFEL8P7PHEUI8DeBpIJKFblPzMoJdvSAzx9EaavcQYeq8OvT2+1S3qbTCaB9rMwFG/tkf/vgEts57EN6iw9Dj/FtQePwZIPJ0nCcd2DUNkgi3e/+MMec5FsCFED936tjMflpD7fK6ZTt6cEabfsi1xbXOsWXLFuzYsQMVp52G0bOn4ca+BXiPhuCgYrAonZKt3EwGc7P3zxhLDccCOMs8eTmejmBDBMSWCNDqwenNtU6rWY+5K75BWAh4iVB+dDG27A5q9sTDQsTdKOzevRsPPPAA/vrXv2LAgAFYv349fD4fnp51p6V53lRXk3MzGczN3j9jLDUMNzNx5KREFwJ4DEAvAA0A6oQQ6vs2KnSmQi5anAhCZop8yGKLeegVJ3l11TYs+6o+7hiXlZd2VE7TUuL34Z0byvGXv/wFs2bNwr59+zBlyhTMmDEDJSXxr9fourhZRMUNne31MpZNzBZycSWAJ6qzB3CnPpS1djlTE1u1TOu5sdniSl4iPDxhqO5NAwF4aHgTJkyYgLFjx2LmzJkYPDh6ow5ldTZC9LaisdfFaCc3t/ZZd1I2vibGOgOzAZyH0DOIU4lJZudF1eZvNbPXNYI3cGhefeZFQ6JqiwshEPzqY7QF92Hg6WMxfvy5WL58OU455ZS4Y8TWZ4+9DY29Lnpzwtmasc1LwRjLbmZ2I2NJsHNbRqcSk7TmRYsLJMOyn4nOqcoB8uEJQ+GTvDgY2IjvX67CrgX3YX/tP7G/OYT+dyzCzUsPxF2zmtpAVPDWorwueuVB9W6MGGMsXXEAd5Dd2zI6VaNaa7vKMScekdBzzZAD5JBuzSj+6DF899KtCO0J4Mgx1+PIy6vR2Nyqec306rMrKa+L1mvUKxfLGduMsXTGAdxBdvfs9IJQsvIV+137fRLGDy/BgtWBqJuPytfWYtiMt6NGEyrKSjB+eEnHnt1WauvtaAhi586dWL/ifcyYMQON332D0p9ciFaKfo2xe1Wbma+PvS56G4nw5h2MsUzEc+AOsrtn58SyJLXEuIOtbXhr7bdxNx+hsOiY25Z7xqu21mPB6sCheWyD87UdbMLelf9AW0sQQy76PU499VRs374dRUVFAMzNVRuJ3QtcpjUnzJt3MMYyEQdwBzmxFtfuxCStUQIzS8qCoXDHGm8jItyK/Wv/hYZlr6CtqQGFx52B7XsOKEqoRgK43jXT2xpVrs9+f8UQw7bE4s07GGOZiAO4gxLp2aV66Y/Z5WNazATvg4GNqP/nI2ip34G8o05A8fi7kNd7YMf5lRnfetds6rw6zXP4CySM6NM94dfBGduMsUzDAdxBVnt2Ti5nir0xOGtQLyzduCupYwKRNd1aQbwtdBAeKQ/eLsVAbgH6Tb4X4ZKyuF3olEu+9K6ZvOZbzZ6mUFYs/WKMMbO4kEsaMSo2kigrldasuqy8FAtWB6KO3fbDVjS+/wJaWltx+IR7IRBZ4623fWxshTc1Zl5HsteKMcbcxoVcMpBTy5n05o6TUSB5cH/FEIzo0x3Vizdh27ZtOPjxK/hhzdvo1q0bDiwSYrQAAA9jSURBVDt1AtpEG4g8hnu/m8kLUPbOtXrivPSLMdZZ8DKyNOLUcianglow1AYgEljvOakVP8z+HRrXL8XNN9+Mx2o+AIZe0LHFpx4CTGd8V5SVYFnVKJTw0i/GWCfHATyNOLXO20pQs7KO+/BCL7788ksAQHl5OaZMmYLPP/8cP5k8FQ/+Z7vp800uL7U8b+3kmnjGGMsEPAeeZmpqA5jx5gbsaYqst/b7JJw/9Ags3bgr4cx0s3Pgfp+E6eMGR51fjWgLI7TpfbR+/Ap6FnfD+vXr4fUeCqZ6m6NIXkJhbg4ag6Gks+x5sw7GWDbiOfAM1tw+NA1ENgV5afm2jq/NZKarBbaZFw3RnTsGgOnjBndkgqtlrS/5bCe+ql2G/f99AU3fbcbw4cPx0EMPRQVvQH/IvkteDhqakg/eAC/9Yox1bhzA04yZhDO9Hci0lqLNvGhIR3b2tJr1UTcFQKRnrKQWHN/KeQtj77gH/fr1w4Nz52LChAnweA7NwshBX2tMh4COnn227PjFGGNu4QCeImaHe80mnMmlRWOPOePNDYZbjqqt/w6FhepNwebNm/HFF19g9OjROOecc/Dcc89h8uTJyM3NjXpt09/YoLuFaOx+3WrtUh6Ph8YZY0wfJ7GlgJVdycwmnPkLpLhjVr66VnPuWnljYGYJ1q5du3DjjTdi0KBBuPbaaxEOh5GTk4MpU6bEBe/bX1+vG7xL/D7NXnnsDYvdO7gxxli24gCeAlZ2JTOzPadP8kIIxG820qadkCjfGNTUBjQzzXv7fThw4AAeeOAB9O/fH48//jh+9atf4cMPP4yb55aZGfK3suyL9+ZmjDFzeAjdokSGd40KtMQec/zwkqisc7nsqfKcenXB1Rw42Ip+VQvhIVLtDctrsVesWIFp06bhggsuwMyZM3Hcccfpvn4ztdT7Vi2E3ydB8hJC4UNnV1v2xXtzM8aYORzALUi0VrneDltqx5yzfJvhzlpGGeWx5CHu2LrlQggEv1iO1obvcAsRJp1yFNatW4chQ+LPrdZWtbltrfNLHkJxgaSbhW73Dm48n84Yy1YcwC3QG97VCwp6O2ypHVMAmLN8G0b06a553MrRA1H56lrdYXMjzds/Q8O7s3Ew8CmkXn3RdcQ4vLR8G15aDgDbUFwg4Z6xg6NKmKq11axQm0BBbg5q7z6743tqy9Via6snWqDFyc1hGGPMbTwHbkGiw7sVZSWYedEQlPh9IESSumZeNAQVZSWazxWA6rxvTW0AI2ctwdR5dWjVKMLjNag7Hmr4Djv/8QC+n1OJ1oZv0X309TjiV/8L8kTPc+9pCqHytbUdCWR2DGMrj6GWsLZgdQDjh5eoXiureD6dMZbNuAduQTLDu1pFR/TmkbUytDuCkkb3t00I/aHtcCsObluPotMvQ7cRFfDk5mu2W7m8zOyctx7ltdIKsEs37rJlRzGeT2eMZTPugVvgRP3tytEDdbPClczuKtbb74t6btvBJjS8/yJ+ePNPAACpx5Eo+Z/n4T/tUt3gLZPXnDe1tBo+FoBmFn3stXI6wDq1OQxjjKUDDuAW6A2FJ3PMyeWlcUHcSoa22vMqRw9EvqcNe1e/icBTv0HjR/MgRBtEOBKEPVI+CiRzb7+85lyvPrpMvibysjF5OF/tWjkdYHnDE8ZYNuMhdIucqL+t3FNbL1taawjbS4Q2IaKet2bNGjS8+Hvs2b4VeaUnovjMKcg7YkDHc3ySFw9eNASrttZjzopt0NvTRm3NORBfXU0OjmavkV5ynx2UyXechc4Yyza8G1kGmVazHnOWb4sLmsqe7d69e9GtWzfs3r0bP/3FGOwdOAae0jKQIrFN3nVMGcjU6qMDwGXlpXHnVCrx+7CjIQh/gQQhYHmXMV7mxRhj0Xg3sixTUxvAgtWBqEBKAMYPj/R2161bh9tuuw07d+7EypUr0aNHD/gnPIgDKj32wrycuCAprzmfu+IbhIWAt31N+P0VQ7B04y7Vnn+J34dlVaOSWq7FO4oxxlhieA7cZvIyr35VCzFy1hLbanhrrcH+1/JPcOWVV2LYsGFYsWIFJk2ahHA48jgrSWLyDYJc6CUsBBasDqCmNmA4l8zLtRhjLPW4B24jJwuHqAXd5u2fYvUrd+KTHA9uvfVW3H777SguLu74uZVlb3pBWF7SpTXUzcu1GGMs9TiA28hMpbZE53zlYNwWOojWPTuQ+6N+yDtiAHqfdiFm3H4rZq/dh7I/fggvEcJCoMRiVTOjIKw31G13+VPGGGPGeAjdRmY2LUl0q8ybf34MDn66BDue+S12vnoPRGsIBfn5mHLzNDz84e6OACoPgVutapbMki5ersUYY6nHPXAbGfVEE6mlLoTAokWLMKOqCt+tX48uRw5E4elX4sie3TRrqSuPbbaqWTJLuni5FmOMpR4HcBsZBcFE5oo/+OADjBkzBv3798e8efNwySWXRC0JM9pW1Ow8dLJBmLPJGWMstTiA28goCJqdK/7yyy+xZs0aTJgwAT/5yU8wb948VFRUIDc3V/W5evXJrcxDcxBmjLHMwQHcZnpB0KiHvnPnTtx777146qmnUFxcjHHjxiE/Px+5A0birEc+UL0pUDum2rG1cCEVxhjLTJzElkJatdR/PqAI9957L/r3748nn3wSv/nNb7Bu3Trk5+cbJr4pjwno1x6PlUxSHWOMMXdxKdU0sHbtWpSVleHCCy/Egw8+iIEDD/WaR85aolsFLRlOHpsxxlhi0rqUKhFVAxgLoAXAVwCmCCEa3GiLG4QQ+Mc//oE1a9bg/vvvx9ChQ/H555/jmGOOiXus3UVSlEPmWrduXICFMcbSn1tD6O8AOEEIcSKAzwHc7lI7Uu6DDz7AyJEjMX78eNTU1CAYjARLteAN2LvlZuyQuRYuwMIYY+nPlQAuhHhbCNHa/uVyAEe60Y5U2rp1Ky644AKcfvrp2Lp1K5555hnU1dXB59MPlnYWSdFbM57ssRljjKVWOmShXwVgntYPiegaANcAQGlpaaraZLvc3FysWrUKDz74IG688UYUFBSYep6dRVL0hsYJ4Cx0xhjLII4lsRHRvwEcrvKjO4UQ/6/9MXcCGAHgImGiIZmexNbS0qK6ljtVOGmNMcbSn+tJbEKIn+v9nIiuBHA+gJ+ZCd7ZwIngbWUddzLlUhljjKUXt7LQzwFwG4AzhBBNbrQhG1jdvpRrljPGWPZwZR04EX0JIA/A7vZvLRdC/NboeZk+hG43HhJnjLHs4/oQuh4hhPqaKWaJ3WvEGWOMZQ4upZrB7FwjzhhjLLNwAM9gdq4RZ4wxllk4gGeQmtoARs5agn5VCzFy1hIAUN0chZPSGGMs+6VDIRdmglbG+cyLhnDCGmOMdULcA88QamVQg6EwqhdvcqlFjDHG3MQBPENwxjljjDElDuAZgjPOGWOMKXEAzxCccc4YY0yJk9gyBJdBZYwxpsQBPINUlJVwwGaMMQaAh9AZY4yxjMQ9cJtZ2d6TMcYYSxQHcBtZ3d6TMcYYSxQPoduIi60wxhhLFQ7gNuJiK4wxxlKFA7iNuNgKY4yxVOEAbiMutsIYYyxVOInNRlxshTHGWKpwALcZF1thjDGWCjyEzhhjjGUgDuCMMcZYBuIAzhhjjGUgDuCMMcZYBuIAzhhjjGUgDuCMMcZYBiIhhNttMI2IdgHY6nY7ktATwA9uNyKL8fV1Fl9fZ/H1dVYmXd8+QoheRg/KqACe6YholRBihNvtyFZ8fZ3F19dZfH2dlY3Xl4fQGWOMsQzEAZwxxhjLQBzAU+tptxuQ5fj6Oouvr7P4+jor664vz4EzxhhjGYh74IwxxlgG4gDOGGOMZSAO4ClGRNVEtJGI1hHRP4jI73absgkRXUJEG4iojYiyasmIW4joHCLaRERfElGV2+3JNkT0HBHtJKJP3G5LNiKio4hoKRF91v7ZcKPbbbILB/DUewfACUKIEwF8DuB2l9uTbT4BcBGA991uSDYgIi+AvwI4F8DxACYR0fHutirrPA/gHLcbkcVaAdwihDgOQDmA67Lld5gDeIoJId4WQrS2f7kcwJFutifbCCE+E0JscrsdWeRkAF8KITYLIVoAvALgApfblFWEEO8DqHe7HdlKCPGtEGJN+7/3AfgMQIm7rbIHB3B3XQVgkduNYExHCYBvFF9vR5Z8+LHOh4j6AigDsMLdltgjx+0GZCMi+jeAw1V+dKcQ4v+1P+ZORIZ25qSybdnAzPVltiGV7/HaU5ZxiKgLgAUAbhJC7HW7PXbgAO4AIcTP9X5ORFcCOB/AzwQvxLfM6PoyW20HcJTi6yMB7HCpLYwlhIgkRIL3HCHE6263xy48hJ5iRHQOgNsAjBNCNLndHsYMrAQwgIj6EVEugEsBvOFymxgzjYgIwLMAPhNCPOJ2e+zEATz1HgfQFcA7RFRHRE+63aBsQkQXEtF2AKcCWEhEi91uUyZrT7i8HsBiRJJ/5gshNrjbquxCRHMBfARgIBFtJ6Jfu92mLDMSwOUARrV/5tYR0XluN8oOXEqVMcYYy0DcA2eMMcYyEAdwxhhjLANxAGeMMcYyEAdwxhhjLANxAGeMMcYyEBdyYSxDEVEPAP9p//JwAGEAu9q/Prm9dnmq27QYwMXtNacZYw7iZWSMZQEimg5gvxDiTzHfJ0T+n7c5fP6UnIcxdggPoTOWZYjoGCL6pL1I0BoARxFRg+LnlxLR39r/fRgRvU5Eq4joYyIqVzneb9r3rl/cvi/4NI3zHNFeiMTf/vMp7fveryWi2WbPxxgzh4fQGctOxwOYIoT4LRHp/T//C4CHhBDL23dqegvACSqPO7n9+y0AVhLRWwD2K88DAJGOOEBEQxEpGXyaEKKeiLpbPB9jzAAHcMay01dCiJUmHvdzREp4yl8XE5FPCBGMedxiIcQeACCiGgA/AfAvnfOMAjBPCFEPAPLfFs7HGDPAAZyx7HRA8e82RG8Lmq/4N8Fcwltssoz89YHYByqOq5ZgY/Z8jDEDPAfOWJZrTyzbQ0QDiMgD4ELFj/8N4Dr5CyIapnGYs4nIT0QFAC4AsMzgtP8GcKk8dK4YQjd7PsaYAQ7gjHUOtyEy5P0fRPb4ll0HYGR7stmnAK7WeP4HAF4GUAtgrhCiTu9kQoh1AB4C8D4R1QGotng+xpgBXkbGGNNFRL8BcIIQ4ia328IYO4R74IwxxlgG4h44Y4wxloG4B84YY4xlIA7gjDHGWAbiAM4YY4xlIA7gjDHGWAbiAM4YY4xloP8PJ7vJdqIghukAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7, 4))\n",
    "plt.scatter(y_train, y_train_pred_lr)\n",
    "plt.plot([-2, 2], [-2, 2], '--k')   #数据已经标准化，3倍标准差即可\n",
    "plt.axis('tight')\n",
    "plt.xlabel('True price')\n",
    "plt.ylabel('Predicted price')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "#岭回归／L2正则\n",
    "#class sklearn.linear_model.RidgeCV(alphas=(0.1, 1.0, 10.0), fit_intercept=True, \n",
    "#                                  normalize=False, scoring=None, cv=None, gcv_mode=None, \n",
    "#                                  store_cv_values=False)\n",
    "from sklearn.linear_model import  RidgeCV\n",
    "\n",
    "#设置超参数（正则参数）范围\n",
    "alphas = [ 0.01, 0.1, 1, 10,100]\n",
    "#n_alphas = 20\n",
    "#alphas = np.logspace(-5,2,n_alphas)\n",
    "\n",
    "#生成一个RidgeCV实例\n",
    "ridge = RidgeCV(alphas=alphas, store_cv_values=True)  \n",
    "\n",
    "#模型训练\n",
    "ridge.fit(X_train, y_train)    \n",
    "\n",
    "#预测\n",
    "y_test_pred_ridge = ridge.predict(X_test)\n",
    "y_train_pred_ridge = ridge.predict(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE: 0.181042313407\n",
      "RMSE: 0.425490673701\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print (\"MSE:\",metrics.mean_squared_error(y_test, y_test_pred_ridge))\n",
    "print (\"RMSE:\",np.sqrt(metrics.mean_squared_error(y_test, y_test_pred_ridge)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEKCAYAAAA4t9PUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VfWd//HXh4SwhDUQtoRVNkH2iNqqiIJ7XVlsa6utjt38Tau/mam/sdPO2N9sddZ2utkZp9PaGUlALSquiLW1orlhk0VWxXsJEDaBBCEk+cwf94TepjfhxuRuyfv5ePDw3HO+95zPPULeOdvnmrsjIiLyUXVJdwEiIpLdFCQiItImChIREWkTBYmIiLSJgkRERNpEQSIiIm2iIBERkTZRkIiISJsoSEREpE1y011AKgwcONBHjRqV7jJERLJKRUXFQXcvPNu4ThEko0aNIhQKpbsMEZGsYma7ExmnU1siItImChIREWkTBYmIiLSJgkRERNpEQSIiIm2iIBERkTZRkIiISJsoSEREOqDIkRP8/fPvUHX8ZNK3pSAREemAykIRfvSrnZyu96RvS0EiItLBNDQ4SysiXDx2IEX9eiR9ewoSEZEO5vWdB9nzwYcsKhmeku0pSEREOpjSUIR+Pbty5eTBKdleUoPEzK42s61mtsPMHoiz/H4z22xmG8xspZmNDOZPN7M3zGxTsGxxzHvMzP7azLaZ2RYz++NkfgYRkWzywYlaXti0j5umF9EtNycl20xa918zywG+D8wHIkC5mS13980xw9YCJe5+wsy+BHwHWAycAD7r7tvNbBhQYWYvuPsHwJ3AcGCiuzeY2aBkfQYRkWzzy3WV1NY1pOy0FiT3iGQ2sMPdd7l7LfA4cGPsAHdf5e4ngpergeJg/jZ33x5MVwJVQGNP/C8BD7l7Q7C8KomfQUQkqywpD3NeUR8mDeuTsm0mM0iKgHDM60gwrzl3Ac81nWlms4E8YGcw6xxgsZmFzOw5MxvXTvWKiGS1jXuOsnnvMRan8GgEkhskFmde3Buazex2oAR4uMn8ocDPgc81HoEA3YCT7l4C/AR4tJl13hOETejAgQMf8SOIiGSP0lCYvNwu3DCtpd/Z218ygyRC9FpGo2KgsukgM5sHPAjc4O6nYub3AZ4FvuHuq5usd1kw/SQwNd7G3f0Rdy9x95LCwrN+U6SISFY7ebqep9bu4ZrzhtC3Z9eUbjuZQVIOjDOz0WaWB9wGLI8dYGYzgB8TDZGqmPl5REPiZ+5e1mS9TwGXB9NzgG1Jql9EJGu8sGkfx07WpfQie6Ok3bXl7nVmdi/wApADPOrum8zsISDk7suJnsrqBZSZGcD77n4DsAi4FBhgZncGq7zT3dcBfwf8wszuA6qBu5P1GUREskVpKExx/x5cNGZAyredtCABcPcVwIom874ZMz2vmfc9BjzWzLIPgOvasUwRkawWPnyC13cc4r554+nSJd7l6eTSk+0iIlmurCKCGSwoKU7L9hUkIiJZrL7BWRoKc8m4wpQ0aIxHQSIiksVe33GQyqMnWZSmoxFQkIiIZLXSUJh+Pbsyf1JqGjTGoyAREclSR2pqeXHT/pQ2aIxHQSIikqV+uW4PtfWpbdAYj4JERCQLuTtLQhGmFPVNaYPGeBQkIiJZaOOeY2zZe4xF56f3aAQUJCIiWak0FKZbbhdumDYs3aUoSEREss3J0/U8tS5o0NgjtQ0a41GQiIhkmRc27eN4mho0xqMgERHJMkvKwwwv6MGFaWjQGI+CREQki4QPn+C3Ow+xcNbwtDRojEdBIiKSRcpC4WiDxlnpa4nSlIJERCRL1Dc4SysiXDqukGFpatAYj4JERCRL/OZMg8bMuMjeSEEiIpIlSkNh+vfsyrxJg9Jdyu9RkIiIZIEjNbW8tGk/N81Ib4PGeBQkIiJZ4KmgQePiDGiJ0pSCREQkw7k7S8rDTC3uy8Qh6W3QGI+CREQkw7295yjv7DuecRfZGylIREQyXGODxk9kQIPGeBQkIiIZ7OTpen65rpJrpwzNiAaN8ShIREQy2PMbow0aF5ZkzpPsTSlIREQy2JLyMCMKenLh6Mxo0BiPgkREJEO9f+gEb+w6xKKS4oxp0BiPgkREJEOVVYTpYnBrBjVojEdBIiKSgc40aBxfyNC+mdOgMR4FiYhIBvr19gPszcAGjfEoSEREMlBZKEJBfh7zzh2c7lLOSkEiIpJhDtfU8uLmfdw0vYi83Mz/MZ35FYqIdDJPrt3D6XrPyAaN8ShIREQyiLtTFgozrbgvE4b0Tnc5CVGQiIhkkA2RoEFjlhyNgIJERCSjlIbCdO+auQ0a41GQiIhkiA9r61m+rpJrzxtKn+6Z2aAxHgWJiEiGeH7TXo6fqmNhFjw7EktBIiKSIZaUhxk5oCcXjilIdymtoiAREckAuw/VsHrXYRaVDMcscxs0xpPUIDGzq81sq5ntMLMH4iy/38w2m9kGM1tpZiOD+dPN7A0z2xQsWxznvd8zs+pk1i8ikiploUi0QePMzG7QGE/SgsTMcoDvA9cAk4BPmtmkJsPWAiXuPhVYCnwnmH8C+Ky7TwauBv7FzPrFrLsE6IeISAfQ2KBxzvhChvTtnu5yWi2ZRySzgR3uvsvda4HHgRtjB7j7Knc/EbxcDRQH87e5+/ZguhKoAgrhTEA9DPxZEmsXEUmZ17YfYN+x7GjQGE8yg6QICMe8jgTzmnMX8FzTmWY2G8gDdgaz7gWWu/veljZuZveYWcjMQgcOHGhV4SIiqVRaHqYgP48rsqBBYzzJDJJ4V4s87kCz24ESokcasfOHAj8HPufuDWY2DFgIfO9sG3f3R9y9xN1LCgsLW128iEgqHKo+xctb9nPzjOxo0BhPbhLXHQFij9OKgcqmg8xsHvAgMMfdT8XM7wM8C3zD3VcHs2cAY4EdwV0NPc1sh7uPTc5HEBFJrsYGjdl6WguSGyTlwDgzGw3sAW4DPhU7wMxmAD8Grnb3qpj5ecCTwM/cvaxxvrs/CwyJGVetEBGRbOXulIbCTBveL2saNMaTtOMod68jej3jBWALUOrum8zsITO7IRj2MNALKDOzdWa2PJi/CLgUuDOYv87MpierVhGRdFgfOcq2/dUszuKjEUjuEQnuvgJY0WTeN2Om5zXzvseAxxJYf6+21igiki6NDRqvnzY03aW0SXZe2RERyXIf1tbz9LpKrp2SXQ0a41GQiIikwXMbow0as/kieyMFiYhIGiwpDzNqQE8uGJ1dDRrjUZCIiKTYewdrePPdwyzMwgaN8ShIRERSrKwinLUNGuNRkIiIpFBjg8bLJgzKygaN8ShIRERS6LVtB9h/7BSLSjrG0QgoSEREUmpJeZgB+XlcPjE7GzTGoyAREUmRjtCgMZ6O80lERDLck2v3UNfgLDo/+58diaUgERFJAXdnSXmY6cP7MX5w9jZojEdBIiKSAuvCH7C9qprFHexoBBQkIiIpURqK0KNrDtdPze4GjfEoSEREkuxEbR1Pr482aOyd5Q0a41GQiIgk2XNv76P6VF2HPK0FChIRkaRbEgozemA+54/qn+5SkkJBIiKSRO8erOGtdw+zsKS4QzRojEdBIiKSRGWhjtWgMR4FiYhIktTVN7C0IsLcCYMY3KdjNGiMR0EiIpIkr20/QNXxUyzsAN+C2BIFiYhIkiwpDzOwVx5XnDso3aUklYJERCQJDlafYuWWKm6eUUTXnI79o7ZjfzoRkTR5ck3QoLGDn9YCBYmISLtzd0pDYWaM6Me4DtagMR4FiYhIO1vb2KCxExyNgIJERKTdlYXC9Oiaw3UdsEFjPAoSEZF2FG3QuJfrpnbMBo3xKEhERNrRig7eoDEeBYmISDsqLQ8zZmA+JSM7ZoPGeBQkIiLtZNeBat567zALS4Z32AaN8SQcJGZ2sZl9LpguNLPRyStLRCT7lFVEyOli3DqzKN2lpFRCQWJm3wK+Dvy/YFZX4LFkFSUikm3q6htYVhFh7oRCBnXgBo3xJHpEcjNwA1AD4O6VQMd/ykZEJEG/2tY5GjTGk2iQ1Lq7Aw5gZvnJK0lEJPs0Nmi8fGLHbtAYT6JBUmpmPwb6mdkfAS8DP0leWSIi2ePA8VO88k4Vt8ws7vANGuPJTWSQu/+Dmc0HjgETgG+6+0tJrUxEJEs8uTYSNGjsuN+C2JKEgiQ4lfWKu79kZhOACWbW1d1PJ7c8EZHMFm3QGGHmiH6MHdQ5Lx0negz2GtDNzIqIntb6HPDTZBUlIpIt1rz/ATuqqjvVk+xNJRok5u4ngFuA77n7zcCks77J7Goz22pmO8zsgTjL7zezzWa2wcxWmtnIYP50M3vDzDYFyxbHvOcXwTo3mtmjZtY5mtmISEYqC4XpmZfDdVOHpbuUtEk4SMzsIuDTwLPBvBZPi5lZDvB94BqiofNJM2saPmuBEnefCiwFvhPMPwF81t0nA1cD/2Jm/YJlvwAmAlOAHsDdCX4GEZF2VXOqjqfXV3LdlKH06pbQlYIOKdEg+SrwAPCEu28Knmp/5SzvmQ3scPdd7l4LPA7cGDvA3VcFRzoAq4HiYP42d98eTFcCVUBh8HqFB4C3Gt8jIpJqz769l5ra+k59WgsSvNhO9AihgehRxe2AETxT0oIiIBzzOgJc0ML4u4Dnms40s9lAHrCzyfyuwGeIhpyISMqVhcKMKcxnVidq0BhPokHyC+BPgI1EAyUR8TqWxQ2fIJxKgDlN5g8Ffg7c4e5Nt/sD4DV3/3Uz67wHuAdgxIgRCZYsIpKYnQeqKX/vCA9cM7FTNWiMJ9EgOeDuT7dy3REg9nivGKhsOsjM5gEPAnPc/VTM/D5Er8d8w91XN3nPt4ie6vpCcxt390eARwBKSkrOdvQkItIqZaFog8ZbOlmDxngSDZJvmdm/AyuBMz/s3f2JFt5TDowLrqfsAW4DPhU7wMxmAD8Grnb3qpj5ecCTwM/cvazJe+4GrgKuiHOUIiKSdHX1DSxbE2HuhEEM6t25GjTGk2iQfI7onVJd+d2pLQeaDRJ3rzOze4EXgBzg0eBC/UNAyN2XAw8DvYCy4NDwfXe/AVgEXAoMMLM7g1Xe6e7rgB8Bu4E3gvc84e4PJfg5RETa7NWtBzhw/FSnfZK9qUSDZJq7T2ntyt19BbCiybxvxkzPa+Z9j9FMm3p377z32IlIRlgSCjOwVzfmdsIGjfEkevvv6jjPgIiIdDpVx0/yyjtV3DqzqFM2aIwn0d/uLwbuMLN3iV4jMcCDBwlFRDqNJ9fsob7BO+X3jjQn0SC5OqlViIhkgWiDxjCzRvZn7KBe6S4nYyTaRn53sgsREcl0a94/ws4DNXzn1nPSXUpG0Qk+EZEElZZHggaNQ9NdSkZRkIiIJKDmVB3PbKjk+qlDye/EDRrjUZCIiCTg2Q1q0NgcBYmISAJKgwaNM0d07gaN8ShIRETOYkdVNaHdR1hcMrzTN2iMR0EiInIWZRVhcroYN6tBY1wKEhGRFpyub2BZxR4un6gGjc1RkIiItODVrQc4WH2KRXqSvVkKEhGRFiwpD1PYuxtzJxSmu5SMpSAREWlG1fGTrNpaxS0zi8hVg8Zmac+IiDTjiaBBo05rtUxBIiISR2ODxpKR/TmnUA0aW6IgERGJo2L3EXYdqGGRnmQ/KwWJiEgcS8rD5OflcN0UNWg8GwWJiEgT1afqePbtvVw/dZgaNCZAQSIi0sSzGyo5UVuv01oJUpCIiDRRGopwTmE+M0f0S3cpWUFBIiISY0fVcSp2H2Hx+WrQmCgFiYhIjLJQhNwuxs0zitNdStZQkIiIBE7XN7BsTYTLJw6isHe3dJeTNRQkIiKBVe9UcbC6Vk+yt5KCREQkUBqKNmi8TA0aW0VBIiICVB07yaqtB7h1ZrEaNLaS9paICLDsTINGXWRvLQWJiHR67k5ZKMzsUQWMUYPGVlOQiEinF9p9hF0Ha1ioo5GPRE1kWvDEmgh7j55kQH4eBfl5DOiVx4D8bhT0yqN3t1w9rCTSQZxp0DhVDRo/CgVJC57ZsJdX3qmKuywvpwsFvxcweRTkd4uZVvCIZIPqU3U8u2EvN04fRs88/Uj8KLTXWvDonedz8nQ9h2pqOVxdy8GaUxyuruVwTdPpWt47VMPh6lpqauvjrqul4CnIj76OTkfnKXhEUuOZ9ZV8eFoNGttCQXIW3bvmUNSvB0X9eiQ0PpHgOZRg8PTP78qA3wsbBY9IeysNhRk7qBczhqtB40elIGlnbQmeQzWnOBQTNoeqTyUcPF1zLE7Y/G66ID+PgUHwFOTn0ae7gkdkR9Vx1rz/AQ9ee67+PbSBgiTNPkrwHK6p5VAQPL+bruVwEESHamrZfegEh6pPnTV4CvK7BQHTGDbdFDzSaZQ2NmicWZTuUrKagiTLdO+aw7B+PRiWhOA5XFNL9am6uOuJDZ7fnVZrPMXW7cz0yAH5anYnWeF0fQNPrIlwxbmDGNhLf2fbQkHSwX3U4DlcU8vB4NRadDoaPI3T7x9uPnimDe/HlZMGc+WkwYwd1EtHMpKRXlGDxnaT1CAxs6uBfwVygH93979rsvx+4G6gDjgAfN7dd5vZdOCHQB+gHvhrd18SvGc08DhQAKwBPuPutcn8HJ1JW4NnU+UxXty0j4df2MrDL2xl9MB85gehMmNEf3K6KFQkM5SWhxnUuxtzxqtBY1uZuydnxWY5wDZgPhAByoFPuvvmmDFzgTfd/YSZfQm4zN0Xm9l4wN19u5kNAyqAc939AzMrBZ5w98fN7EfAenf/YUu1lJSUeCgUSsrnlPj2HT3JS1v289Lm/byx8yCn650B+Xlcce4grpw0hIvHDaR715x0lymd1P5jJ7nob1fyhTnn8PWrJ6a7nIxlZhXuXnK2cck8IpkN7HD3XUFBjwM3AmeCxN1XxYxfDdwezN8WM6bSzKqAQjM7ClwOfCpY/F/AXxI9epEMMqRvdz5z4Ug+c+FIjp08za+2HuDFzft57u19lIYi9Oiaw6XjBzJ/0hCumDiI/vl56S5ZOpFlayI0ODqt1U6SGSRFQDjmdQS4oIXxdwHPNZ1pZrOBPGAnMAD4wN0bT8xHgu1IBuvTvSufmDaMT0wbRm1dA6t3HeKlzdGjlRc27aeLwfmjCrhy8hCunDSY4QU9012ydGDRBo0RZo8uYPTA/HSX0yEkM0jinQyPex7NzG4HSoA5TeYPBX4O3OHuDRb/qm1z67wHuAdgxIgRrShbkikvtwuXji/k0vGFPHTjZN7ec5SXNu/nxU37+fYzm/n2M5uZOKQ3V04azPxJQzivqI8u1ku7Kn/vCO8erOErc8emu5QOI5lBEgFijxuLgcqmg8xsHvAgMMfdT8XM7wM8C3zD3VcHsw8C/cwsNzgqibtOAHd/BHgEotdI2v5xpL2ZGVOL+zG1uB//98oJ7D5UEw2Vzfv5t1U7+O4rOxjatzvzJw1m/qTBXDB6AHm5algtbbOkPEyvbrlcO2VIukvpMJIZJOXAuOAuqz3Abfzu2gYAZjYD+DFwtbtXxczPA54EfubuZY3z3d3NbBWwgOidW3cAv0ziZ5AUGjkgn7svGcPdl4zhcE0tK4OL9aWhMD97Yze9u+cyd8Igrpw8mDnjC+ndvWu6S5Ysc/zkaVa8vZebZqhBY3tK2p509zozuxd4gejtv4+6+yYzewgIufty4GGgF1AWnL54391vABYBlwIDzOzOYJV3uvs64OvA42b2/4G1wH8k6zNI+hTk57GwZDgLS4bzYW09v9lxkBc37WPlO1UsX19J1xzjonMGBqfABjO4T/d0lyxZ4JkNe6MNGnWRvV0l7fbfTKLbfzuO+ganYvcRXtq8jxc372f3oROAHoKUxNz8g9epPlnHi/ddqr8jCUj09l8FiWQtd2d7VXVwsX4f6yNHARg1oCdXTh7C/EmDmamHICWwff9x5v/za3zjunO5+5Ix6S4nK2TCcyQiSWVmjB/cm/GDe/OVuWN/7yHI/3z9XR55bdeZhyDnTxrCJXoIslMrDYXJ7WLcNENPDLQ3BYl0GPEegnypyUOQl4wbyJWTh3D5xEEU6CHITqO2roEn1uxh3rmD1aAxCRQk0iE1fQjyzXcP8eKm/WduL+5iUDKqILiuMoQRA/QQZEf2yjtVHKqpZdH5xekupUPSNRLpVNz9zEOQL23ezzv7jgMwcUjvoLmkHoLsiD7/03I2VR7l9a9fTm6OnkVKlK6RiMTR9CHI9w+d4MXgDrDvr9rB9/QQZIez7+hJXt1axRfnnKMQSRIFiXRqIwb0TOghyPmTBnPZBD0EmY3UoDH5FCQiAT0E2fFEGzSGuWB0AaPUoDFpFCQicfTIyzlzeqvpQ5DfeGoj33hqI9OK+555XmWcHoLMSG+9e5j3Dp3g/1w+Lt2ldGi62C7SCi09BDl/0mCunDxED0FmkPtL1/Hipv2UPziPHnl6hqi1dLFdJAlaegjyp799j5/8+l09BJkhGhs03jyjWCGSZAoSkTY420OQ3bt24dJxhcyfNJgrzh2shyBT6On1ezl5uoHF5+sie7IpSETaiR6CzCyloTDjB/diWnHfdJfS4ekaiUiSNfcQ5ITBvbnqvCHcOrOIkQN0R1F72rb/OFeqQWOb6RqJSIZo6SHI772yne+u3M7s0QUsmFXMdVOGkt9N/yzbqrQ8TNcc42Y1aEwJHZGIpNHeox/yxJo9LK2I8O7BGnrm5XDtlKEsmFXMBaMLdEvxR1Bb18CFf7uSC0YX8MPbZ6W7nKymIxKRLDC0bw++MncsX77sHCp2H6EsFOGZDZUsrYgwoqAnC2YVc8vMIor763pKolZu2c/hmlo9yZ5COiIRyTAnaut4fuM+ykIR3th1CDP42DkDWDhrOFdNHqJbWc/ic//5Flv2Huf1By7X8zxtpCMSkSzVMy+XW2YWc8vMYsKHT7BsTYSlFRG+tmQdvbvlcv20oSyYNZyZI/rp1FcT+46e5FfbDvDly8YqRFJIQSKSwYYX9ORr88bzx5eP4813D1NWEeaptZX8z1thxhTmR099zShmSF/1/YLfNWhcWKLvHUklndoSyTLVp+pYsWEvZRVhyt87QheDS8YVsrCkmHnnDu60T9I3NDhz//FVhvbtzuP3XJTucjoEndoS6aB6dctl0fnDWXT+cN47WMPSigjL1kS497/X0rdHV26YNoyFJcVMKerbqU59vfXeYXYfOsFXr1CDxlRTkIhksVED8/mTqyZw3/zx/HbnQcpCEUpDYX6+ejcTBvdmwaxibppRRGHvjv895aXlYXp3y+Wa84amu5ROR0Ei0gHkdDEuGVfIJeMKOfrhaZ7ZUElZKMJfr9jC3z3/DnMnFLJg1nAunzioQ37j47GTp1mxcS+3zFSDxnRQkIh0MH17dOXTF4zk0xeMZEfVccoqIjyxZg8vb6miID+PG6cPY+Gs4Uwa1ifdpbabp9dXRhs06tmRtNDFdpFOoK6+gV9vP0hZRZiXN1dRW9/A5GF9WDCrmBunF2V9V+Ibv/86J2vref5rl3Sq60LJpovtInJGbk4X5k4cxNyJgzhSU8vy9ZWUVYT5q6c38zcrtjDv3MEsmFXMnPGF5OZk16mvrfuOsz78AX9x/SSFSJooSEQ6mf75edzxsVHc8bFRbNl7jKUVEZ5au4fnNu6jsHc3bplRxIJZxYwb3DvdpSZkiRo0pp1ObYkItXUNrNpaxdKKCKveqaKuwZk2vB8LZhVzw9Rh9O3ZNd0lxlVb18AFf/MyF50zgB98Wg0a25tObYlIwvJyu3DV5CFcNXkIB6tP8dTaaEfiv3hqI99+ZjNXTR7CglnFXDx2YEa1Hnl5y36OnDjNQl1kTysFiYj8noG9unH3JWO46+LRbKo8RlkozC/XV/L0+kqG9u3OLTOLuHVmMWMKe6W7VEpDYYb27c6l4wrTXUqnpiARkbjMjPOK+nJeUV/+/LpzWbmlirJQmB++upPvr9pJycj+LCwp5topQ+ndPfWnvvYe/ZDXth3gK3PVoDHdFCQiclbdcqNfuHXtlKHsP3aSJ9fuoSwU5uvL3uYvl2/mmvOGsKCkmAtHD6BLin6oL6sIGjTO0mmtdFOQiEirDO7TnS/OOYcvXDqGteEPWFoR4el1lTyxdg/F/Xtw68xiFswqZnhB8r6Mq6HBKQ1FuGjMAEYM0Jd+pZuCREQ+EjNj5oj+zBzRn29eP4kXNu1jaUWE776ynX9duZ0LxxSwcNZwrpkyhJ557fuj5s13D/P+4RPcN18NGjOBbv8VkXa154MPeTL4Mq73Dp0gPy+H66YOZWHJcEpG9m+XhwbvW7KOl7fsp/zBeZ22bX4q6PZfEUmLon49uPfycXxl7ljK3zvC0oowz27YS2kowqgBjd9DX8ywfj0+0vqPnTzNirf3smBWsUIkQyhIRCQpzIzZowuYPbqAb31iMs9t3MfSijD/8OI2/vGlbVw8diALZhVz1eQhrQqE5esqOVXXwOLzdZE9UyhIRCTp8rvlsmBW9CL8+4dOsHRNhGUVEb76+Dp6d8/lhmnDWDCrmOnDz/499GWhMBOH9GZKUd8UVS9nk9TubGZ2tZltNbMdZvZAnOX3m9lmM9tgZivNbGTMsufN7AMze6bJe64wszVmts7MfmNmY5P5GUSkfY0Y0JP754/n1382l/+++wLmnTuYZWsi3PyD3zL/n1/jR7/aSdWxk3Hf+86+Y6yPHGVRyXA1aMwgSbvYbmY5wDZgPhAByoFPuvvmmDFzgTfd/YSZfQm4zN0XB8uuAHoCX3D362Pesw240d23mNmXgdnufmdLtehiu0hmO37yNM9u2MvSigih3UfI6WLMGV/IglnFXHHuILrlRk99/dXTm3hs9W7e/PN5Wd/6PhtkwsX22cAOd98VFPQ4cCNwJkjcfVXM+NXA7THLVprZZXHW60DjN/L0BSrbt2wRSbXe3bty2+wR3DZ7BLsOVLNsTYRlFXv48jtr6NezKzdNL+LG6cN4au0erpw0RCGSYZIZJEVAOOZ1BLighfF3Ac8lsN67gRVm9iFwDLgw3iAzuwe4B2DEiBGJ1CsiGWBMYS/+9KqJ3D9/Ar/ZcZClFRH++633+elv3wNgYUlxeguUP5DMIIl3AjPueTQzux3O9JYuAAAIoElEQVQoAeYksN77gGvd/U0z+1Pgn4iGy+9vyP0R4BGIntpKtGgRyQyNp7fmjC/k6InTLN9QyZ4jH3KJGjRmnGQGSQSIvT+vmDinocxsHvAgMMfdT7W0QjMrBKa5+5vBrCXA8+1Trohkqr49u/KZC0eefaCkRTLv2ioHxpnZaDPLA24DlscOMLMZwI+BG9y9KoF1HgH6mtn44PV8YEs71iwiIq2UtCMSd68zs3uBF4Ac4FF332RmDwEhd18OPAz0AsqCW/ned/cbAMzs18BEoJeZRYC73P0FM/sjYJmZNRANls8n6zOIiMjZqdeWiIjElejtv0l9IFFERDo+BYmIiLSJgkRERNpEQSIiIm2iIBERkTbpFHdtmdkBYPdHfPtA4GA7ltNeVFfrqK7WUV2t01HrGunuZ20l0CmCpC3MLJTI7W+pprpaR3W1jupqnc5el05tiYhImyhIRESkTRQkZ/dIugtohupqHdXVOqqrdTp1XbpGIiIibaIjEhERaRMFSRNm9rCZvWNmG8zsSTPr18y4q81sq5ntMLMHUlDXQjPbZGYNZtbsXRhm9p6ZvW1m68ws6Z0qW1FXqvdXgZm9ZGbbg//2b2ZcfbCv1pnZ8nhj2qmeFj+/mXUzsyXB8jfNbFSyamllXXea2YGYffQHXyKXpLoeNbMqM9vYzHIzs+8GdW8ws5kZUNNlZnY0Zl99M9k1BdsdbmarzGxL8G/xq3HGJHd/ubv+xPwBrgRyg+m/B/4+zpgcYCcwBsgD1gOTklzXucAE4FWgpIVx7wEDU7i/zlpXmvbXd4AHgukH4v1/DJZVp2AfnfXzA18GfhRM3wYsyZC67gT+LVV/n2K2eykwE9jYzPJriX41txH9uu03M6Cmy4Bn0rCvhgIzg+newLY4/x+Tur90RNKEu7/o7nXBy9VEv9mxqdnADnff5e61wOPAjUmua4u7b03mNj6KBOtK+f4K1v9fwfR/ATcleXstSeTzx9a7FLjCgi/pSXNdaeHurwGHWxhyI/Azj1oN9DOzoWmuKS3cfa+7rwmmjxP9sr+iJsOSur8UJC37PNEUb6oICMe8jvCH/+PSxYEXzazCzO5JdzGBdOyvwe6+F6L/0IBBzYzrbmYhM1ttZskKm0Q+/5kxwS8yR4EBSaqnNXUB3BqcDllqZsPjLE+HTP03eJGZrTez58xscqo3HpwSnQG82WRRUvdXMr+zPWOZ2cvAkDiLHnT3XwZjHgTqgF/EW0WceW2+/S2RuhLwcXevNLNBwEtm9k7wm1Q660r5/mrFakYE+2sM8IqZve3uO9taWxOJfP6k7KOzSGSbTwP/4+6nzOyLRI+aLk9yXYlIx/46mzVEW4pUm9m1wFPAuFRt3Mx6AcuAr7n7saaL47yl3fZXpwwSd5/X0nIzuwO4HrjCgxOMTUSA2N/MioHKZNeV4Doqg/9WmdmTRE9ftClI2qGulO8vM9tvZkPdfW9wCF/VzDoa99cuM3uV6G9z7R0kiXz+xjERM8sF+pL80yhnrcvdD8W8/AnR64aZICl/p9oi9oe3u68wsx+Y2UB3T3oPLjPrSjREfuHuT8QZktT9pVNbTZjZ1cDXgRvc/UQzw8qBcWY22szyiF4cTdodP4kys3wz6904TfTGgbh3mKRYOvbXcuCOYPoO4A+OnMysv5l1C6YHAh8HNiehlkQ+f2y9C4BXmvklJqV1NTmPfgPR8++ZYDnw2eBupAuBo42nMtPFzIY0Xtcys9lEf74eavld7bJdA/4D2OLu/9TMsOTur1TfYZDpf4AdRM8lrgv+NN5JMwxYETPuWqJ3R+wkeoon2XXdTPS3ilPAfuCFpnURvftmffBnU6bUlab9NQBYCWwP/lsQzC8B/j2Y/hjwdrC/3gbuSmI9f/D5gYeI/sIC0B0oC/7+vQWMSfY+SrCuvw3+Lq0HVgETU1TX/wB7gdPB36+7gC8CXwyWG/D9oO63aeFOxhTWdG/MvloNfCxF++pioqepNsT83Lo2lftLT7aLiEib6NSWiIi0iYJERETaREEiIiJtoiAREZE2UZCIiEibKEhEWmBm1W18/9LgqfmWxrxqLXROTnRMk/GFZvZ8ouNF2kJBIpIkQa+lHHffleptu/sBYK+ZfTzV25bOR0EikoDgieCHzWyjRb/vZXEwv0vQCmOTmT1jZivMbEHwtk8T80S9mf0waBC5ycz+qpntVJvZP5rZGjNbaWaFMYsXmtlbZrbNzC4Jxo8ys18H49eY2cdixj8V1CCSVAoSkcTcAkwHpgHzgIeD9iG3AKOAKcDdwEUx7/k4UBHz+kF3LwGmAnPMbGqc7eQDa9x9JvAr4Fsxy3LdfTbwtZj5VcD8YPxi4Lsx40PAJa3/qCKt0ymbNop8BBcT7YJbD+w3s18B5wfzy9y9AdhnZqti3jMUOBDzelHQ2j83WDaJaFuLWA3AkmD6MSC2AV/jdAXR8ALoCvybmU0H6oHxMeOriLaqEUkqBYlIYpr7kqmWvnzqQ6I9tDCz0cCfAOe7+xEz+2njsrOI7WF0KvhvPb/7t3sf0R5n04ieYTgZM757UINIUunUlkhiXgMWm1lOcN3iUqLNFX9D9IufupjZYKJft9poCzA2mO4D1ABHg3HXNLOdLkS7/wJ8Klh/S/oCe4Mjos8Q/frcRuPJjO7P0sHpiEQkMU8Svf6xnuhRwp+5+z4zWwZcQfQH9jai30x3NHjPs0SD5WV3X29ma4l2h90FvN7MdmqAyWZWEaxn8Vnq+gGwzMwWEu3OWxOzbG5Qg0hSqfuvSBuZWS+PfiveAKJHKR8PQqYH0R/uHw+urSSyrmp379VOdb0G3OjuR9pjfSLN0RGJSNs9Y2b9gDzg2+6+D8DdPzSzbxH9buz3U1lQcPrtnxQikgo6IhERkTbRxXYREWkTBYmIiLSJgkRERNpEQSIiIm2iIBERkTZRkIiISJv8L5TEVN3D8AwSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha is: 10.0\n"
     ]
    }
   ],
   "source": [
    "mse_mean = np.mean(ridge.cv_values_, axis = 0)\n",
    "plt.plot(np.log10(alphas), mse_mean.reshape(len(alphas),1)) \n",
    "\n",
    "#这是为了标出最佳参数的位置，不是必须\n",
    "#plt.plot(np.log10(ridge.alpha_)*np.ones(3), [0.28, 0.29, 0.30])\n",
    "\n",
    "plt.xlabel('log(alpha)')\n",
    "plt.ylabel('mse')\n",
    "plt.show()\n",
    "\n",
    "print ('alpha is:', ridge.alpha_)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:1094: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().\n",
      "  y = column_or_1d(y, warn=True)\n",
      "/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:491: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.\n",
      "  ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:491: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.\n",
      "  ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:491: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.\n",
      "  ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:491: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.\n",
      "  ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:491: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.\n",
      "  ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:491: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.\n",
      "  ConvergenceWarning)\n",
      "/anaconda3/lib/python3.6/site-packages/sklearn/linear_model/coordinate_descent.py:491: ConvergenceWarning: Objective did not converge. You might want to increase the number of iterations. Fitting data with very small alpha may cause precision problems.\n",
      "  ConvergenceWarning)\n"
     ]
    }
   ],
   "source": [
    "#### Lasso／L1正则\n",
    "#class sklearn.linear_model.LassoCV(eps=0.001, n_alphas=100, alphas=None, fit_intercept=True, \n",
    " #                                   normalize=False, precompute=’auto’, max_iter=1000, \n",
    "  #                                  tol=0.0001, copy_X=True, cv=None, verbose=False, n_jobs=1,\n",
    "   #                                 positive=False, random_state=None, selection=’cyclic’)\n",
    "from sklearn.linear_model import LassoCV\n",
    "\n",
    "#设置超参数搜索范围\n",
    "#alphas = [ 0.01, 0.1, 1, 10,100]\n",
    "\n",
    "#生成一个LassoCV实例\n",
    "#lasso = LassoCV(alphas=alphas)  \n",
    "lasso = LassoCV()  \n",
    "\n",
    "#训练（内含CV）\n",
    "lasso.fit(X_train, y_train)  \n",
    "\n",
    "#测试\n",
    "y_test_pred_lasso = lasso.predict(X_test)\n",
    "y_train_pred_lasso = lasso.predict(X_train)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE: 0.181551727641\n",
      "RMSE: 0.426088872938\n"
     ]
    }
   ],
   "source": [
    "print (\"MSE:\",metrics.mean_squared_error(y_test, y_test_pred_lasso))\n",
    "\n",
    "print (\"RMSE:\",np.sqrt(metrics.mean_squared_error(y_test, y_test_pred_lasso)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从以上实验的三个方法中可以看出岭回归模型的RMSE的值最小，说明岭回归模型对于这个例子效果最好。岭回归的目标函数是在OLS上加了一个正则项."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
